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Abstract. 

We have developed a new empirically-based transport algorithm for use in our 
GSFC two-dimensional transport and chemistry model. The new algorithm contains 
planetary wave statistics, and parameterizations to account for the effects due to 
gravity waves and equatorial Kelvin waves. As such, this scheme utilizes significantly 
more information compared to our previous algorithm which was based only on zonal 
mean temperatures and heating rates. The new model transport captures much of 
the qualitative structure and seasonal variability observed in long lived tracers, such 
as: isolation of the tropics and the southern hemisphere winter polar vortex; the well 
mixed surf-zone region of the winter sub-tropics and mid-latitudes; the latitudinal and 
seasonal variations of total ozone; and the seasonal variations of mesospheric H 2 0. 
The model also indicates a double peaked structure in methane associated with the 
semiannual oscillation in the tropical upper stratosphere. This feature is similar in 
phase but is significantly weaker in amplitude compared to the observations. The model 
simulations of carbon- 14 and strontium-90 are in good agreement with observations, 
both in simulating the peak in mixing ratio at 20-25 km, and the decrease with altitude 
in mixing ratio above 25 km. We also find mostly good agreement between modeled 
and observed age of air determined from SF 6 outside of the northern hemisphere polar 
vortex. However, observations inside the vortex reveal significantly older air compared 
to the model. This is consistent with the model deficiencies in simulating CH 4 in the 
northern hemisphere winter high latitudes and illustrates the limitations of the current 
climatological zonal mean model formulation. The propagation of seasonal signals in 
water vapor and C0 2 in the lower stratosphere showed general agreement in phase, 
and the model qualitatively captured the observed amplitude decrease in C0 2 from the 
tropics to midlatitudes. However, the simulated seasonal amplitudes were attenuated 
too rapidly with altitude in the tropics. Overall, the simulations with the new transport 
formulation are in substantially better agreement with observations compared with our 
previous model transport. 
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1. Introduction 

Two-dimensional (2D) chemistry and transport models have been used extensively 
to study various natural physical processes in the middle atmosphere, and to assess 
the long-term impact of man-made halogens on stratospheric ozone [WMO, 1991, 
1995]. Of fundamental importance in 2D modeling is the proper representation of the 
zonally averaged transport fields which are comprised of mean advective and eddy flux 
components. Previous theoretical development has shown that the residual circulation of 
the transformed Eulerian mean (TEM) formulation provides a reasonable approximation 
to the net advective mass transport in the meridional plane [e.g., Andrews and McIntyre, 
1976; Dunkerton, 1978]. Earlier models which used a TEM or isentropic formulation 
for the mean transport and a constant diffusion coefficient to parameterize the eddy 
flux effects were able to reproduce some general features of the observed zonal mean 
temperature, zonal wind, and long lived tracer fields [e.g., Garcia and Solomon, 1983; 
Ko et al., 1985]. Other model studies have specified various eddy diffusion rates to 
investigate the exchange of mass between the tropics and mid-latitudes [e.g., Weisenstein 
et al., 1996]. 

Subsequent studies noted that planetary wave processes simultaneously generate 
both a mean circulation and eddy flux effects so that a model formulation should account 
for these components in a self-consistent manner. For 2-D models with empirically based 
’fixed’ transport schemes, these quantities were typically diagnosed from empirically 
determined temperatures and net heating rates [e.g., Newman et al., 1988; Yang et al., 
1991; Patten et al., 1994; Fleming et al., 1995]. Here, the wave drag and diffusion were 
not derived explicitly, but were inferred from the residual needed to balance the zonal 
momentum equation. Fully coupled models, in which the zonal mean dynamical fields 
are computed interactively with the ozone and constituent distributions [e.g., Garcia 
and Solomon, 1983], now usually include a parameterization to explicitly compute the 
wave induced drag and diffusion due to planetary wave breaking [Garcia, 1991]. This 
also allows for feedback between the planetary wave field and the zonal mean flow. 

In addition to planetary waves, parameterizations have been developed to compute 
the effects of small scale waves. For example, the parameterization of Dunkerton [1979] 
computed the momentum deposition due to Kelvin wave absorption in the tropical 
upper stratosphere, which is important for simulation of the semiannual oscillation 
(SAO). Lindzen [1981] developed a parameterization to account for the effects of 
breaking gravity waves in the mesosphere. This scheme has subsequently been widely 
incorporated into chemical and dynamical models [e.g., Holton and Zhu, 1984; Garcia 
and Solomon, 1985; Gray and Pyle, 1987; Brasseur et al., 1990]. Other gravity wave 
drag schemes [e.g., Bacmeister, 1993; Pierrehumbert, 1987; Hines, 1997] have recently 
been developed and incorporated into 2D chemical and/or dynamical models [Mengel et 
al., 1996; Summers et al., 1997]. 

These paramaterizations have enabled models to better simulate the general zonal 
mean temperature and zonal wind structure in the middle atmosphere, as well as 
characteristic dynamical features observed in long lived tracers such as the degree of 
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tropical and polar vortex isolation, the sub-tropical/mid-latitude surf zone of the winter 
hemisphere, and the upper stratospheric SAO double-peaked structure [e.g., Gray and 
Pyle, 1987; Brasseur et al., 1990; Choi and Holton, 1991; Yang et al., 1991; Garcia 
et al., 1992; Bacmeister et al., 1995; Nightingale et al., 1996; Jackman et al., 1996; 
Summers et al., 1997]. 

In addition to the development of parameterizations of wave processes, recent 
work in 2D chemical modeling in general has tended to focus on the development 
of fully interactive models [e.g., Garcia et al., 1992; Ko et al., 1993; Bacmeister et 
al., 1995; Summers et al, 1997]. While empirical ’fixed transport’ models provide a 
somewhat simpler framework for studying the atmospheric chemical response to various 
perturbations, interactive models have an advantage in that they can compute the 
feedback between chemistry and transport. They can therefore simulate future chemical, 
radiative, and dynamical changes induced by chlorine loading and other anthropogenic 
perturbations. Interactive models may ultimately be needed to fully understand the 
effects caused by such perturbations. However, they do not at present completely resolve 
all details of the observed zonal mean temperature and zonal wind fields. This point 
also has implications for 2D model parameterizations of eddy transport processes which 
can be quite sensitive to the model representation of the zonal mean flow. Most if not 
all interactive models also have a limitation in that they do not account for the effects 
of synoptic scale waves which are important in and around the tropopause region and 
in the very lower stratosphere. 

In the present study, we have developed an upgraded empirical transport algorithm. 
This utilizes several quantities derived from observations: zonal mean temperature, 
zonal wind, net heating rates, and Eliassen-Palm flux diagnostics for planetary and 
synoptic scale waves. We also account for the effects of gravity waves and equatorial 
Kelvin waves by incorporating the parameterizations mentioned above in which the 
zonal mean flow is constrained to observations. These quantities are then used to derive 
self-consistent residual circulation and diffusion fields. The unique approach here is that 
we have incorporated these wave parameterizations and explicit diagnostic quantities 
into one comprehensive, empirically-based transport scheme. Such a scheme should be 
able to resolve more of the details of middle atmospheric transport compared to our 
previous empirically-based algorithms. 

We have recently used this new transport algorithm in the GSFC 2D chemistry 
model to study seasonal and long-term changes in stratospheric ozone [Jackman 
et al., 1996] and mesospheric water vapor [Chandra et al., 1997], and to study 
stratospheric tracer correlations [Dessler et al., 1998]. We have also used the model 
in several international assessment studies of stratospheric ozone, including the World 
Meteorological Organization (WMO), the Intergovernmental Panel on Climate Change 
(IPCC), and NASA’s Atmospheric Effects of Stratospheric Aircraft Program (AESA). 

In conjunction with these assessments, we have participated in the Models and 
Measurements II intercomparison project (MMII) [Park et al, 1999] in which diagnostic 
tests have been performed separately on the model transport and chemistry components. 

The purpose of this paper is to provide a detailed description and validation of 
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the new transport algorithm. We will show model simulations of standard long-lived 
tracers (CH 4 ,H 2 0, and total ozone) compared with observations. This will help assess 
how well the new model resolves the general zonal mean transport processes in the 
middle atmosphere. More recent emphasis has been put on model assessments of 
potential long term ozone changes due to perturbations in the lower stratosphere. 
Such perturbations include anthropogenic chlorine loading [WMO, 1995]; the injection 
and dispersion of volcanic aerosols [e.g., Solomon et al., 1996; Jackman et al., 1996]; 
and the projected emissions of supersonic aircraft [e.g., Stolarski et al., 1995]. Since 
these ozone simulations are greatly influenced by the model dynamics in the lower 
stratosphere, w T e will address more detailed aspects of the transport in this region. This 
will include model/measurement comparisons of inert radioactive tracers carbon- 14 and 
strontium-90, age of air derived from SF 6 , and the propagation of seasonal cycles in 
water vapor and C0 2 . 

2. GSFC-2D Model 

The 2-D model at NASA/Goddard Space Flight Center (GSFC) was originally 
described in Douglass et al. [1989], and extended to mesospheric heights by Jackman et 
al. [1990]. The model has recently undergone significant development in several areas, 
including improvements in the method of computing the photolytic source term and 
the photolysis of 0 2 [Jackman et al., 1996]. We have also recently updated the model 
reaction rates and photolysis cross sections to the Jet Propulsion Laboratory (JPL) 
1997 recommendations [DeMore et al., 1997]. Improvements to the model transport 
algorithm were summarized in Jackman et al. [1996]. In this paper, we present a 
thorough analysis of the new transport fields and comparisons with observations. A 
detailed description of the transport formulation and methodology are contained in 
the Appendix. We note that the transport has undergone some minor changes since 
the recent MMII intercomparison project, so that the results shown here will differ 
somewhat from those presented in the MMII publications [Park et al., 1999; Hall et al., 
1998]. 

In section 3.2, we present model simulations of H 2 0 compared with observations. 
Our method for computing water vapor was previously described [Fleming et al., 1995], 
and is briefly summarized here. Photochemical production and loss follows the JPL-97 
recommendations. We specify a boundary condition at 660 mb that varies in latitude 
and season based on the water vapor climatology of Oort [1983]. We include a loss 
due to rainout throughout the troposphere which is based on the model climatological 
temperature field. This rainout loss is inversely proportional to the temperature, and is 
largest at the tropical tropopause. There is also a loss due to the sedimentation of ice 
particles formed by heterogeneous chemical processes [Considine et al., 1994]. 
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3. Results 

We will now examine the spatial structure and seasonal evolution of several model 
simulated tracer fields and show comparisons with observations. We will first show 
CH 4 ,H 2 0, and total ozone from our standard model which has been run for 20 years 
to obtain a seasonally repeating steady state simulation. We will then show results 
from time dependent simulations of C0 2 , SFe, and radioactive products of atmospheric 
nuclear bomb tests, carbon-14 and strontium-90, to illustrate more detailed aspects of 
the model transport in the upper troposphere and lower stratosphere. We will also 
compare our new model simulations with those using the previous version of the model 
transport to illustrate how the new algorithm improves the tracer simulations. 

3.1. Methane 

Figure 1 shows latitude-height sections of methane from the model simulation 
(solid) and observations (dashed). Here, we use the combined climatology compiled by 
Randel et al [1998] based on data from the Halogen Occultation Experiment (HALOE) 
and the Cryogenic Limb Array Etalon Spectrometer (CLAES) instruments aboard 
the Upper Atmosphere Research Satellite (UARS). Here, the CLAES data have been 
used to fill in the high latitude regions during winter where HALOE does not obtain 
measurements. This climatology has been averaged using potential vorticity (PV) as 
a horizontal coordinate and then mapped to equivalent latitude (the latitude of an 
equivalent PV distribution arranged symmetrically about the pole - see Randel et al. 
[1998] for details). This averaging procedure increases the effective latitudinal coverage 
of the data. Results are similar to standard latitude averaging except at high northern 
latitudes during winter and spring when the vortex can be significantly displaced from 
the pole. 

The structure and seasonal variations of the model compare qualitatively well with 
the data. Characteristic dynamical features are evident, such as the isolation of the 
tropics and the well mixed surf-zone region of the winter subtropics and midlatitudes. 
The fact that the simulated gradients in the tropics and subtropics compare reasonably 
well with the data supports the use of the HRDI wind data to derive the model wave 
driving and diffusion quantities in this region (section A. 2. a). The comparison in Figure 
1 also shows some discrepancies. For example, the absolute values of CH 4 in the model 
tend to be underestimated in the middle and upper stratosphere. A similar difference 
was noted in comparisons between the HALOE+CLAES climatology and the MACCM2 
general circulation model [Waugh et al., 1997]. At high latitudes, the model captures 
the isolation of the southern hemisphere (SH) vortex in July and October in the lower 
and middle stratosphere. By contrast in the upper stratosphere, flat tracer gradients 
illustrating strong horizontal mixing into the polar region are evident during October in 
both the data and model (the 300 ppbv contour). Note also that high latitude descent 
during autumn in the upper stratosphere is evident in both hemispheres. Small CH 4 
values characteristic of mesospheric air extend down to 35-40 km in both the model and 
data. 
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The model shows some differences at high northern latitudes, which can be 
more clearly seen in the month-height cross sections in Figure 2. Here, the northern 
hemisphere (NH) has been shifted by 6 months to facilitate the visual comparison of the 
polar regions during winter and spring. In both hemispheres, descent within the vortex 
is indicated by the steady decline of the observed CH 4 isopleths from autumn through 
early spring, followed by a sharp increase indicating the spring break-up of the polar 
vortex and mixing in of midlatitude air. However, the model does not fully resolve the 
observed low values of CH 4 within the northern polar vortex. This is not surprising since 
this region experiences large longitudinal variability. This varibility is accounted for in 
the equivalent latitude mapping of the data, but is not resolved with the standard zonal 
mean model formulation. In contrast, the southern polar region typically exhibits less 
variability, so that the model agrees fairly well with the data in this region in Figures 1 
and 2. 

Figures 3 and 4 provide a heuristic illustration of the processes that control the long 
lived tracer distributions in the stratosphere. Figure 3 shows the standard model CH 4 
distribution in color for January and September, overlaid with streamlines depicting 
the residual circulation. Upward motion is present in the tropics throughout the year, 
providing the advective flux of mass from the troposphere to the stratosphere. Strong 
descent is evident in the winter polar region of both hemispheres, with ascent in the 
summer upper stratosphere. The summer lower stratosphere in both hemispheres is 
characterized by very weak vertical motion. Figure 4 shows the same CH 4 distributions 
in the solid contours, with the model K yy values overlaid in color. Weak mixing is 
evident in the tropics and summer hemisphere, except for the very lower summer 
stratosphere where decaying baroclinic wave activity from the troposphere generates 
mixing below 20 km. Strong mixing is seen in the surf zone region of both hemispheres 
(20° to 50°), which is reflected in the flat gradients of the CH 4 field. In the winter 
polar regions, the mixing is stronger in the NH compared with the SH, reflecting the 
hemispheric asymmetries in planetary wave activity. Note also the minimum in mixing 
at 60°S in September extending from ~20 to 40 km. This coincides with the maximum 
in u at the vortex edge and results in sharp horizontal tracer gradients with much larger 
values of methane outside compared to inside the vortex. 

Referring back to Figure 1, the UARS data for April shows the well-known 
’double-peaked’ structure at low latitudes in the upper stratosphere. This feature is 
associated with the meridional circulation induced by the SAO and is most pronounced 
during the NH spring. This characteristic is also seen in the month-latitude sections of 
the data and model CH 4 at 2.2 mb (Figure 5). The model, which utilizes the Kelvin 
wave parameterization discussed in section A.2.c, captures the overall time evolution 
and phasing seen in the observations, but has weaker gradients compared to the data. 
Similar findings were reported in other modeling studies of the SAO using climatological 
data [Choi and Holton, 1991], and results from the NCAR MACCM2 [Randel et al., 
1994; Waugh et al., 1997]. Choi and Holton [1991] discuss the possible limitations of 
simulating the strength of the SAO circulation using climatological data as is done here. 
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3.2. Water vapor 

Figure 6 shows latitude-height sections of H 2 O from the model simulation (solid) 
and UARS observations (dashed). In the stratosphere and lower mesosphere, we use 
the combined climatology compiled by Randel et al. [1998]. This is based on data from 
HALOE augmented with data from the Microwave Limb Sounder (MLS) in the polar 
regions. As was done for CH 4 , the H 2 O climatology is mapped in equivalent latitude 
coordinates. For the mesosphere above 60 km, Figure 6 includes HALOE version 18 
H 2 0 data for October 1991 through December 1996 averaged in standard latitude 
coordinates. Again the model reproduces the overall structure and seasonal variations 
seen in the data. Mixing ratios in lower tropical stratosphere are ~3-4 ppmv, consistent 
with the flux of dry air into the stratosphere from the tropical tropopause. The oxidation 
of methane causes the water vapor concentrations to increase with height above the 
lower stratosphere, to maximum of 6-7 ppmv in the summer upper stratosphere and 
lower mesosphere. 

The seasonal variation of H 2 0 in the equatorial lower stratosphere is shown 
in Figure 7. Here we show the HALOE data for 1993-1995, along with the model 
climatological simulation repeated for every year with no interannual variations. The 
model qualitatively captures the annual cycle at the tropical tropopause. Dry (moist) 
air enters the stratosphere during the NH winter (summer), coincident with the 
coldest (warmest) tropopause temperatures. This seasonal variation is advected slowly 
upwards by the residual circulation, and maintains its structure for a long time period 
because of the weak horizontal and vertical diffusion characteristic of the tropical lower 
stratosphere [e.g., Mote et al., 1996]. In the upper stratosphere above 35 km, the 
model indicates a semiannual cycle which is generally consistent with the HALOE data. 
The model amplitudes are underestimated relative to HALOE, although the signature 
in the data is contaminated by interannual variability not represented in the model. 
This semiannual feature is related to the SAO circulation which brings down water 
vapor-rich air from the lower mesosphere during the equinoxes. The model circulation 
and the accompanying signature in the tracers tend to be a bit more apparent during 
March-April than September-October, which is consistent with observations indicating 
that the SAO circulation is stronger during the NH spring [e.g., Holton and Choi, 1988] 
(see also Figures 1 and 5). 

The lower stratospheric seasonal signal is more readily seen in the quantity, 

2CH 4 + H 2 0 (defined as H as in Mote et al. [1998]), which is quasi-conserved since the 
chemical destruction of methane produces two molecules of water vapor. Figure 8 again 
shows that the model qualitatively reproduces the characteristics of the data, with the 
signal extending up to above 10 mb. The amplitude and phase lag based on Fourier 
analysis of the annual harmonic of this signal are shown in Figure 9. Here, the quantities 
plotted are relative to those at the tropical tropopause. The model phase lag variation 
with height shows good agreement with the data, with an average phase propagation of 
about .31 mm/sec between 100 and 10 mb. However, the model amplitude is attenuated 
faster with increasing altitude than observed. This discrepancy was characteristic of the 
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various models analyzed in the MMII intercomparisons [Hall et al., 1998]. 

To understand the transport processes that affect the seasonal signal propagation in 
this region, we have plotted our model transport fields in the lower tropical stratosphere 
in Figure 10. These can be compared with those derived from HALOE data by Mote et 
al. [1998]. Here we present values for the different seasons and the annual average. The 
annual mean residual vertical velocity (w*) is in good agreement with the Mote et al. 
analysis, with a minimum of .22 mm/sec at 19-21 km, increasing with height to about 
.4 mm/sec at 32 km. The average phase propagation of .31 mm/sec in the model and 
data for H in Figure 9 generally reflects the average vertical velocity between 16 and 
32 km. The minimum w* at 19-25 km contributes to a slight increase in the vertical 
gradient of the model phase (i.e., slower propagation) at these altitudes in Figure 9, 
although this feature is weaker than observed. Consistent with previous findings, the 
seasonal variations in the lower stratosphere reveal a maximum w* during NH winter 
(January), and a minimum during the SH winter (July). This is attributed to the 
greater extratropical planetary wave driving during the NH winter compared to the 
SH winter which in turn controls the seasonal variations of the vertical mass flux in 
the lower tropical stratosphere [e.g., Rosenlof, 1995; Mote et al., 1996]. The seasonal 
variation in vertical velocity is thought to be responsible for the seasonal variation in 
temperature in this region, which in turn drives the seasonal signal in water vapor. 

The model vertical diffusion rates in Figure 10 also are generally similar to Mote et 
al., 1998 and the observational analysis of Hall and Waugh [1997b]: values are .01-. 02 
m 2 s ~ 1 at 17-29 km, and increase to greater than .1 m 2 s -1 above 30 km (see Appendix). 

To estimate the dilution rate of midlatitude air into the tropics, we computed 
the annual mean K yy value averaged over 35°S to 35°N, divided by the equator-35 0 
distance squared. The results in Figure 10 reveal a minimum annual mean dilution rate 
(maximum time scale) near 22 km. This is qualitatively consistent with observational 
analyses that suggest that the tropics are considerably more isolated from the 
midlatitudes above about 19 km than below [e.g., Schoeberl et al., 1997; Mote et al., 
1998]. Our maximum time scale near 22 km is about 2 years, a bit longer than reported 
by other previous studies [e.g., Schoeberl et al., 1997; Hall and Waugh 1997b; Randel et 
al., 1998], but much shorter than the 6-7 year time scale derived by Mote et al. [1998]. 
The seasonal variation in model dilution rates in Figure 10 is consistent with variations 
in extratropical wave activity. Maximum dilution rates occur during the NH winter 
(January) and SH late winter-spring (October), with minimum dilution in April and 
July. 

As discussed by Mote et al. [1998] and Hall et al. [1998], entrainment of midlatitude 
air is largely responsible for the attenuation of the observed taperecorder signal in the 
lower tropical stratosphere, although vertical diffusion is also important despite the 
relatively small values. The observed entrainment is a minimum at 19-25 km which 
coincides with both the reduced attenuation of the taperecorder signal and the sharp 
reduction in the vertical gradients of CH 4 and H observed in this region. Our model 
simulation noticeably attenuates the seasonal signal faster with increasing altitude than 
is observed, especially in the 19-25 km region, and shows only a very weak reduction of 
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the vertical gradients of CH 4 and H. 

Reasons for the model overattenuation of the seasonal signal are unclear. The 
HALOE data has been noted to underestimate the H 2 O annual cycle near the tropopause 
[e.g., Mote et al, 1996] so that the signal attenuation with height in Figure 9 may 
be underestimated by HALOE. The model dilution time scales above 20 km are 
quantitatively consistent with some observational studies, but are significantly shorter 
than the Mote et al. [1998] analysis. To test if the model entrainment rates are too large, 
we set K yy to a very small value of lO^m 2 ^" 1 throughout the year at 35°S to 35°N in the 
lower stratosphere. This corresponds to a dilution time scale of 49 years and implies an 
extremely isolated tropical region. However, the corresponding amplitude attenuation 
of the seasonal signal in H was only slightly diminished, and was still significantly 
greater than observed. This suggests that our model entrainment of midlatitude air in 
Figure 10 is reasonable. Increasing the model vertical diffusion may effectively reduce 
the signal attenuation to be closer to observations, however, this will diminish the 
model-measurement agreement in other areas [Hall et al, 1998], For example, the phase 
of the signal will be propagated too quickly and the vertical gradient of the age of air 
will be decreased (see section 3. 5. a, Figure 22). Our model vertical velocity and vertical 
diffusion rates shown in Figure 10 are in fairly good agreement with the observations 
so that model discrepancies in these quantities do not seem to significantly contribute 
to the overattenuation. Numerical diffusion introduced by the advection scheme could 
artificially damp out the signal, however such diffusion is small in the current scheme 
[Lin and Rood, 1996]. It is possible that the overattenuation is caused by the relatively 
coarse vertical resolution of the model (2 km) which may not enable the simulation to 
properly resolve the propagation of seasonal signals. This possibility will be addressed 
in future work. 

3. 2. a. Mesospheric Water vapor 

Returning to Figure 6, both the model and data show the overall sharp decrease 
in H 2 0 with height in the upper mesosphere due to the photo-dissociation by Lyman 
alpha flux. As a result, seasonal variations in the gravity wave-driven mesospheric 
meridional circulation are reflected in the water vapor distribution. In the summer 
extratropics, ascent advects water vapor-rich air upwards from the stratosphere to 
the upper mesosphere, and wintertime descent brings down dry air from the lower 
thermosphere. Net photochemical effects are small in winter, and although there is 
some net photochemical loss of H 2 0 in the summer upper mesosphere, this effect is not 
large enough to offset the water vapor increase due to vertical advection. Transport due 
to vertical diffusion has a secondary effect as disscussed below, and horizontal diffusive 
transport is small in our model simulation. 

The seasonal variations of the model and HALOE water vapor at 80 km and 65 
km for 3 latitude regions are shown in Figures 11 and 12. At both altitudes, the model 
captures the overall seasonal variability observed by HALOE. There has been debate over 
the relative importance of vertical advection and diffusion in controlling the water vapor 
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distribution in the mesosphere [e.g., Garcia and Solomon, 1985; Holton and Schoeberl, 
1988; Smith and Brasseur, 1991; Nedoluha et al., 1996]. A time tendency analysis of our 
model simulations revealed that throughout most of the mesosphere, vertical advection 
has a significantly larger influence on H 2 0 than does diffusion [Chandra et al., 1997]. 
However, diffusion does play a secondary role in determining the extratropical seasonal 
variations seen in Figure 11. Maximum vertical diffusion occurs during the solstices, 
coincident with the strong westerly and easterly midlatitude mesospheric jets (section 
A.2.b). This diffusion acts to increase (decrease) H 2 0 in the upper (lower) mesosphere. 
In the upper mesosphere, the diffusive transport therefore enhances the water vapor 
increase in summer and partially offsets the H 2 0 decrease in winter due to vertical 
advection. This effect maximizes at 80 km, and results in the secondary maxima seen in 
Figure 11 during mid- winter at midlatitudes. This effect is somewhat more pronounced 
in the SH winter because the stronger zonal westerlies, and hence the larger difference 
between the background zonal wind and the gravity wave phase speed (u — c), allow for 
larger diffusion compared to the NH winter (section A.2.b). We also note that these 
secondary maxima are partially generated by the horizontal advection of moist air from 
the summer to the winter hemisphere induced by the meridional circulation. However, 
this effect is substantially less than the contribution due to vertical diffusion. At high 
latitudes, u*, and hence the horizontal advection, are quite small, so that the secondary 
mid-winter maxima are due only to vertical diffusion. 

At the equator in Figure 11, the model at 80 km shows a dominant semiannual 
cycle qualitatively similar to the data, although the absolute magnitudes of the model 
H 2 0 are underestimated by 1-1.5 ppmv. Here, the effect from vertical diffusion is 
small throughout the year, so that the semiannual variation in vertical advection, 
with maximum upward w* during the solstices (Figure A. 7), controls the water vapor 
seasonality. 

The overall seasonal variations in transport become weaker with decreasing altitude 
so that the corresponding seasonality in H 2 0 is much less pronounced at 65 km in 
Figure 12. This is especially true at the equator where a rather flat seasonal structure 
is seen in both the model and HALOE data. 

3.3. Total Ozone 

Our model total ozone simulation and comparison with TOMS version 7 data were 
reported in Jackman et al. [1996]. Since that publication, the simulation has improved 
somewhat in the southern hemisphere midlatitudes and we provide a summary of our 
current model-TOMS comparison here. Figure 13 shows the model total ozone for 1990 
conditions of total chlorine loading, along with the TOMS data averaged over 1988 to 
1992, and the difference (model minus TOMS). Much of the observed seasonal structure 
is qualitatively simulated by the model, including the on-the-pole maximum of ~440 
DU during the NH spring and the minimum in the NH autumn; the seasonality in the 
tropics, with a maximum during October, a secondary maximum during April, and a 
minimum during the NH winter; the off-the-pole maximum during the southern late 
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winter-spring; and the very low ozone at the southern pole characterizing the Antarctic 
ozone hole during spring. 

There are some notable discrepancies in the simulation: (1) the model predicts 
larger amounts of ozone at midlatitudes during summer and early fall, especially in 
the NH; (2) the magnitude of the simulated southern off-the-pole-maximum is not as 
large as observed; and (3) the model does not simulate the near elimination of strong 
latitudinal gradients at high southern latitudes observed following the break-up of the 
ozone hole. The model also shows noticeably less ozone compared to TOMS at high 
southern latitudes throughout most of the year. To check if discrepancy (3) is due 
to an underestimation of horizontal mixing during the spring breakup of the Antartic 
vortex, we performed an additional model run with K yy set to a very large value of 
4xl0 lo cm 2 s _1 poleward of 50°S throughout the lower stratosphere during October and 
November. However, this gave only a marginal improvement compared to the standard 
model simulation. This suggests that other deficiencies in the model are causing the 
larger than observed latitudinal gradients during the austral spring and summer. 

3.4. Carbon-14 and Strontium-90 

The radioactive products of atmospheric nuclear bomb tests, carbon-14 ( 14 C) 
and strontium-90 ( 90 Sr), are very useful in testing the transport in a model [e.g., 
Jackman et ah, 1991; Prather and Remsberg, 1993; Kinnison et al., 1994], 14 C and 90 Sr 
were produced in large quantities in the late 1950s and early 1960s by aboveground 
nuclear tests conducted by the United States and the Soviet Union. These radioactive 
products were measured extensively by aircraft and balloonsondes in the troposphere 
and stratosphere [see Kinnison et al., 1994 and references therein]. The measurements 
after January 1, 1963, are especially important since there were no atmospheric nuclear 
detonations from 1963 through 1966, and so no excess production of 14 C or 90 Sr needs 
to be accounted for in model simulations during this time period. There is a natural 
source of 14 C from galactic cosmic rays which we do not attempt to simulate. All model 
simuations and measurements shown in this section are the excess 14 C above the natural 
background state. 

The distribution of 14 C from 1963 through 1966 is dependent primarily on the 
transport processes of the stratosphere and troposphere as 14 C is in the form of 14 C02 
and acts like a passive tracer. The half-life of 14 C is 5730 years and provides a negligible 
loss over this time period. The distribution of 90 Sr in the middle 1960s is dependent 
on both the transport processes of the stratosphere and troposphere and the settling 
velocity of aerosol particles since 90 Sr rapidly coalesces onto these particles. The half-life 
of 90 Sr is 28 years and provides only a very small loss over this time period. 

The initial conditions for 14 C and 90 Sr input into our model were taken from 
Prather and Remsberg [1993] and their latitude and altitude dependent distributions 
were specified on October 15, 1963, and on October 15, 1964, respectively. Time and 
hemisphere dependent ground boundary conditions for 14 C and latitude and altitude 
dependent settling velocities for 90 Sr were also specified from Prather and Remsberg 



ill '1 



13 


[1993]. 

Our time-dependent simulation of 14 C is shown in Figure 14 for the months of 
January and July in years 1964, 1965, and 1966. The peak at middle to high latitudes 
persists through July 1965 although the model transport moves the enhanced 14 C in 
the tropical lower stratosphere up through the middle (January 1965) and upper (July 
1965) stratosphere before being totally dissipated (January 1966). Strong gradients in 
the Northern lower stratosphere between the tropics and extra-tropics persist through 
January 1966 indicating the very slow tranport and mixing that takes place in this 
region. The final distribution shown in July 1966 is similar to that of a long-lived source 
gas (see CH4 in Figure 1), except the largest amounts of 14 C are found in the upper 
stratosphere with smaller amounts at lower altitudes. It thus takes about 3 years for a 
primarily Northern Hemisphere lower stratospheric source distribution of a substance to 
become relatively well mixed throughout the stratosphere. 

Detailed comparisons between the model and measurements at 31°N are shown 
in Figure 15. Also shown are model results from an older version (designated ”1995 
Model”) as well as the present version of our model transport (designated ”1998 
Model”). The ”1998 Model” is clearly an improvement over the ”1995 Model” during 
the period of the 14 C simulation. The 14 C peak between 20 and 25 km and the very 
large gradient between 10 and 20 km in the measurements are represented much better 
by the ”1998 Model.” 

We also show comparisons of 90 Sr between the two versions of our model and 
the measurements at two latitudes, 9°N and 34°S, in Figures 16 and 17, respectively. 
At both of these latitudes in which the atmospheric transport characteristics are 
substantially different, it is apparent that the ”1998 Model” simulates the 90 Sr more 
realistically than the ”1995 Model.” The two model simulations especially diverge from 
one another in late 1965. The very excellent agreement between the ”1998 Model” and 
the measurements in October 1966 is probably somewhat fortuitous, given that the 
real atmosphere undergoes interannual dynamical variability that is not contained in 
our T998 Model’ transport. It is significant, however, that the modeled 90 Sr generally 
simulates the gradual decrease of measured stratospheric "Sr over 1965 and 1966, while 
representing the peak concentration at 20-25 km in a very reasonable manner. 

There are substantial differences between our ”1998 Model” and ”1995 Model” 
transport algorithms which are reflected in these 14 C and "Sr simulations. While the 
advective transport by the mean circulation exhibits some differences, most of the 
changes seen in Figures 15-17 are due to changes in the meridional and vertical diffusive 
transport. The ”1998 Model” is generally less diffusive in a globally averaged sense. 
However more importantly, the sharp changes in diffusion that occur across certain 
atmospheric regions are much better resolved in the ”1998 Model”. These include the 
changes in vertical mixing across the tropopause, and the changes in horizontal mixing 
between the tropics and subtropics in the stratosphere. 
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3.5. Sulfur hexafluoride 

SFg is a nearly inert anthropogenically-produced tracer. It has had a relatively 
steady source at the ground over the past 2-3 decades, allowing for a quasi-linear 
temporal increase in atmospheric concentration [e.g., Geller et al., 1997]. The time lag 
of SF 6 between a given point in the stratosphere and a reference location at the ground 
determines the mean age of air, T [Elkins et al., 1996; Harnisch et al., 1996]. This 
has become a widely used diagnostic of model transport [e.g., Hall and Plumb, 1994; 
Waugh et al., 1997; Bacmeister et al., 1998; Hall et al., 1998]. The annually averaged T 
derived from SF 6 is very similar to the mean age determined from the first moment of 
the age spectrum in model simulations [Hall and Plumb, 1994; Hall and Waugh, 1997a; 
Hall et al., 1998]. However, seasonal variations of the age determined from SF6 can be 
significant. 

Our model simulation of T determined from SFg for various seasons is shown in 
Figure 18. Here we have taken the area weighted global mean value at the surface as 
the reference point. In the troposphere, the vertical gradients are quite small due to 
the rapid vertical mixing, and there is a 1-1.5 year time lag between midlatitudes of the 
northern and southern hemispheres. This interhemispheric exchange time is in good 
agreement with SFg measurements [e.g., Harnisch et al., 1996; Geller et al., 1997] and 
largely results from the tropospheric horizontal diffusion imposed in the model (see 
section A. 2. a). The very lower stratosphere exhibits sharp vertical gradients in age 
due to the drastic reduction in vertical mixing across the tropopause into the lower 
stratosphere. The age of air contour shapes in the stratosphere are qualitatively similar 
to other long lived tracers (e.g., CH 4 , Figure 1), consistent with the mean circulation. 
The youngest stratospheric air occurs just above the tropical tropopause (about .6 years 
older than the global mean surface value), with the oldest air of ~6 years occurring in 
the mesosphere and polar winter upper stratosphere. The mean age in the mesosphere 
is relatively uniform owing to the vigorous meridional circulation and vertical mixing 
characteristic of this region. The older ages seen in the polar stratosphere in fall and 
winter are indicative of the descent of older mesospheric air within the polar vortex. As 
expected, this feature is more pronounced in the SH where, relative to the NH, the polar 
vortex is more isolated from the younger mid-latitude air. The isolation of the SH polar 
region in the lower stratosphere persists into October, as indicated by the small 5.5 year 
isopleth centered at 50 mb at the south pole in Figure 18d. This is also reflected in the 
CH 4 simulation in Figure Id. 

Observations of T can be made in certain regions using measurements of SF 6 from 
aircraft [Elkins et al., 1996], balloons [Harnisch et al., 1996; Patra et al., 1997; Moore et 
al., 1998], and the space shuttle [Rinsland et al., 1993]. Mean ages can also be derived 
from CO 2 measurements after taking into account the annual cycle in CO 2 at the ground 
and a small correction for CH 4 oxidation [e.g., Boering et al., 1996]. Ages derived from 
CO 2 in this way were found to be similar to those derived from SFg data [Waugh et al., 
1997; Hall et al., 1998]. 

Figure 19 shows the mean age for October/November 1994 (top) and Jan- 
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uary/February 1996 (bottom) computed from three model simulations at 20 km 
compared with that derived from SFg data at 19-21 km taken onboard several NASA 
ER-2 flights [Elkins et al., 1996]. The global mean surface time series of Geller et al. 
[1997] is used as the reference point for the SF 6 measurements. Note that in deriving T 
from both the model and observations, we have not accounted for the small nonlinear 
effect in the time increase of SF 6 . This effect slightly increases the mean age compared 
to the simple time lag method, with the difference approaching .5 years for T of ~6 
years [Volk et al., 1997]. However, for the present analysis we are mainly interested in 
the relative differences in T between the model and measurements. 

The base ‘1998 Model’ simulation (solid line) in Figure 19 captures much of the 
latitudinal variation seen in the data, including the sharp gradient between the equator 
and ±20°. Ages near the equator are ~1.5 years in both the model and data. The 
model-data agreement of the low-latitude gradients supports the use of the HRDI 
wind data to derive the model wave driving and diffusion quantities in the tropics 
(section A. 2. a). While there is a fairly large spread in the observations at middle and 
higher latitudes, the model T tends to be a little younger than observed, especially 
in January/February 1996. For example near 50°N, during January /February 1996, 
the mean of the data is about 5 years whereas the model indicates about 4 years. 
Underestimation of T is a characteristic seen in most stratospheric models [Hall et al., 
1998]. 

Figure 20 shows vertical profiles of T from three different model runs for the various 
latitudes and seasons indicated. These are plotted along with T derived from the balloon 
SF 6 data of Harnisch et al. [1996] (asterisks and squares) and the Observations of the 
Middle Stratosphere campaign (OMS) [Moore et al., 1998] (triangles). Again we have 
used the global mean value at the surface [Geller et al, 1997] as the reference function 
for the data. This gives stratospheric ages that are .25 to 1 year younger compared to 
ages derived using the NH tropospheric reference function used by Harnisch et al. [1996]. 
At low to middle latitudes, the standard T998 model’ (solid line) shows good agreement 
with the observations, including the transition from very weak vertical gradients in the 
troposphere to much stronger gradients in the lower stratosphere. The data generally 
indicate a transition to little or no vertical gradient above about 25 km, and the T998 
model’ qualitatively reproduces this feature in most of the comparisons. At midlatitudes 
above ~22 km, the model agrees well with the data at 44°N but not as well at 35°N. 
These data show age differences of about 1 year in the middle stratosphere, even though 
they are at similar latitudes and seasons (September). The data were taken during 
different years (1996 versus 1993), so some of the difference may be due to interannual 
dynamical variability effects which are not contained in the model. At high latitudes, 
the model agrees well with the late winter observations in March taken outside of the 
polar vortex at 68°N (Figure 20h, squares), and the summer observations at 65°N 
(Figure 20e). The model also shows a region of reduced vertical age gradients at 10-15 
km in the summer profile at 65°N which is somewhat indicated in the data. However, 
the model did not capture the very old air (8-10 years) observed inside the vortex above 
20 km at 68°N in mid- winter (Figure 20g) and late winter (Figure 20h, asterisks). 
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Figures 19-20 also show the mean age from our previous ’1995 Model’ version 
(dotted line) which was discussed in the 14 C and 90 Sr simulations (section 3.4). Relative 
to the observations and the ’1998 model’ simulation, T from the ’1995 model’ is 
significantly underestimated, and the associated latitudinal gradient is much flatter. 
The vertical gradients in the lower stratosphere (Figure 20) are also generally weaker 
than observed. This is consistent with the T995 Model’ simulation of 14 C and 90 Sr 
(section 3.4) and reflects the overly diffusive nature of this previous model transport. 

While these model simulations assume SF6 is an inert tracer, it is possible that 
SF 6 incurs a small photo-chemical loss in the mesosphere. If such a loss exists it would 
effectively increase the model T to be older than the base simulations shown above. 

To test this, we included a loss identical to C2CIF5 (CFC-115) in our ’1998 Model’ 
simulation in Figures 19-20 (dashed line). CFC-115 photolytic losses are taken from 
DeMore et al. [1997] and result in a lifetime of about 600 years in our model. Significant 
increases in the model age occur at middle and high latitudes, and at the higher 
altitudes in the tropics when including this mesospheric loss. At 20 km (Figure 19), T is 
generally .5 to 1 year older at middle to high latitudes, with ages increased by nearly 2 
years at the SH pole during spring. This latter result reflects the mesospheric character 
of air within the SH vortex, even at levels reaching down to the lower stratosphere. 
Including the mesospheric loss has increased Y in the mid-stratosphere by 1 year at low 
latitudes, 2 years at middle latitudes, and 3 years at high latitudes (Figure 20). 

Overall, the model Y is in better agreement with the data at the higher altitudes 
and latitudes when including this mesospheric loss, although the comparison is not as 
good at 44°N (Figure 20f). The comparisons inside the NH vortex are also improved, 
however, the model still underestimates the observations by 1 to 2 years (Figure 20g, 
Figure 20h, asterisks). We note that a loss rate comparable to CFC-115 is considered a 
probable upper limit and is shown here primarily as a qualitative example of how such 
a loss could effect Y. Differences between the base model simulation and observations 
could reflect other limitations of the model such as the grid resolution. Furthermore, 
discrepancies within the NH polar vortex are probably due to the fact that the vortex 
is frequently centered off the pole and exhibits large year to year variability and as 
such cannot be resolved by the current zonal mean climatological model formulation. 

In contrast, the model does a reasonable job in resolving the SH vortex (Figure lc-d, 
Figure 18c-d) which has smaller longitudinal and interannual variability. 

3. 5. a. Age of air sensitivity to K zz 

An important transport component in 2D models is the vertical eddy diffusivity 
{K zz ). In the mesosphere, diffusion coefficients can be determined from gravity wave 
parameterizations [e.g., Lindzen, 1981], which provide for reasonably good long lived 
tracer simulations (section 3. 2. a, see also Garcia et al. [1985]). Determination of K zz in 
the upper troposphere and lower stratosphere is very important for proper treatment 
of the tropopause boundary and troposphere-stratosphere exchange processes. In this 
region of our model, we scale K zz based on the vertical temperature gradient, so that 
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a stronger lapse rate, indicative of more rapid convective overturning, implies a larger 
value of diffusion. This methodology implies a large mixing rate in the troposphere, and 
a very small I\ zz in the lower stratosphere, with a sharp gradient across the tropopause 
boundary (see section A.2.d). We specify a lower limit on K zz of .01-. 02 m 2 s~ 1 in the 
lower stratosphere following the observational analyses of Hall and Waugh [1997b] and 
Mote et al. [1998]. 

An evaluation of the lower stratospheric vertical diffusion specified in models was 
recently made by Hall et al. [1998]. These authors discussed the influence of K zz 
on the propagation of seasonal signals in the tropics and concluded that a value of 
K zz < .1 m 2 s 1 is most realistic. Here, we illustrate the sensitivity of the model global 
distribution of T to the lower stratospheric minimum in K zz . This will help to evaluate 
the qualitative realism of the current model I\ zz field. 

Figures 21-22 show the mean age from the base T998 model’ simulation along 
with the SF 6 observations as in Figures 19-20. We also include two additional model 
scenarios in which the lower stratospheric minimum I\ zz is increased to .1 m 2 s~ l and 
1 ra 2 s _1 , respectively, but are otherwise identical to the base ’1998 model 1 . Increasing 
the minimum rate of vertical mixing reduces the mean age throughout the stratosphere. 
Increasing the minimum to .1 m 2 s~ 1 decreases T only slightly in the tropical lower 
stratosphere, but the difference increases to as much as .5 years in the lower stratosphere 
extratropics and in the middle stratosphere at all latitudes. 

With a minimum of 1 m 2 s ~ l , T is signicantly younger everywhere compared to 
observations and the base model scenario. Maximum differences of 3 years occur at 
middle and high latitudes in the middle stratosphere. In addition, both the vertical 
and horizontal gradients in mean age are drastically reduced in this model run. This 
unrealistic scenario is far outside the range of the observations, and is similar to the 
’diffusive regime’ discussed by Hall et al [1998] in their analysis of seasonal signal 
propagation. Although the model scenario with the minimum K zz set to .1 m 2 s~ 1 is 
within the observational values of T, the overall model-measurement agreement is not 
as good relative to the base model case. This provides evidence that a minimum K zz of 
.01-. 02 m 2 s~ 1 gives the most realistic model simulation, as is suggested by the analyses 
of Hall and Waugh [1997b] and Mote et al [1998]. Our model-measurement comparisons 
of CH4 also reveal that a minimum K zz of .01-. 02 m 2 s~ l gives the best simulation. 

3.6. Carbon dioxide 

The source of carbon dioxide at the ground includes a large secular trend of 
~1.4 ppmv/yr along with a large biosphere-induced seasonal cycle. This seasonal 
variation is largest at high northern latitudes and is strongly latitude dependent. At 
stratospheric levels, CO2 is nearly inert: it has no photochemical loss, and has only a 
small production from CH4 oxidation. The transport characteristics of the stratosphere 
are therefore reflected in the propagation of the CO2 seasonal cycle. As shown in 
previous observational analyses, this signal originates at the tropical tropopause and 
propagates vertically in the tropics with a gradual loss of amplitude, similar to the 
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observed behavior of H . The CO 2 signal is also transported rapidly poleward with only 
a small decrease in amplitude just above the tropopause. At higher altitudes (above 
~460K, ~19 km) the midlatitude signal is very small, indicating that the barrier to 
poleward transport from the tropics becomes increasingly strong with height [Strahan 
et al., 1998; Mote et al., 1998]. Simulation of these characteristics in CO 2 therefore 
provides a diagnostic of model transport. 

Our time dependent model simulation of CO 2 uses lower boundary conditions based 
on monthly mean global surface observations for 1979 to 1995 [Conway et al, 1994]. 
The simulation was continued through 1996 by using the 1995 boundary conditions 
and adding 1.4 ppmv. To account for the small CO 2 source from CH 4 oxidation, the 
simulation includes the corresponding time dependent surface boundary conditions of 
CH 4 [WMO, 1995] and is run with all model photo-chemistry included. 

Figure 23a-b shows the model simulated annual cycle amplitude of CO 2 , defined 
as one-half of the difference from peak to trough of the seasonal variation, and 
the corresponding phase. The amplitude maximizes at the ground in the northern 
hemisphere extratropics (~7.5 ppmv), and decreases equatorward and with increasing 
height in the troposphere. The phase in the troposphere is fairly uniform in the 
tropics and northern hemisphere, reflecting the strong horizontal and vertical mixing 
of the region. In the lower stratosphere, the amplitude attenuates with height and a. 
strong phase shift is seen with increasing latitude and altitude away from the tropical 
tropopause. 

Previous studies have shown that the seasonal cycle in CO 2 in the northern 
hemisphere midlatitude lower stratosphere is propagated from the tropics and not 
from the underlying troposphere [e.g., Strahan et al., 1998]. Figure 24 shows the 
model simulation at 55°N for five pressure levels increasing in altitude from the middle 
troposphere to the lower stratosphere. There is a dramatic decrease in amplitude and 
an increase in phase lag between 282 and 212 mb, illustrating a distinct difference in the 
airmass characteristics of the troposphere and lower stratosphere. This was also seen 
in the sharp gradient in amplitude and phase across the tropopause in Figure 23, and 
suggests that the model reasonably simulates the separation between the troposphere 
and stratosphere. 

Comparisons of the model simulation in certain regions of the stratosphere can 
be made with CO 2 observations from the ER-2 field campaigns [Boering et al, 1996; 
Strahan et al., 1998]]. We first examine the simulation at the tropical tropopause which 
essentially is the boundary condition for CO 2 entering the the stratosphere. Figure 
25 shows time series of the model simulation (solid line) and observations for 1992 
through 1996. Here, the observations are taken as the average of the surface time 
series at Mauna Loa (19°N) and Samoa (14°S) with a two month time lag (dashed 
line). This time series gives an excellent fit to the ER-2 CO 2 data just above the 
tropical tropopause (triangles), and therefore provides a proxy for the stratospheric CO 2 
boundary condition [Boering et al., 1996]. The model agrees fairly well in phase, lagging 
the observations by about 1-2 weeks. However, the model amplitude is underestimated 
by roughly a third (1.2 vs. 1.8 ppmv). While the model captures the observed maxima 
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fairly well, the amplitude discrepancy is due primarily to the model underestimating the 
observed minima. We note that at lower altitudes, the model exhibits a significantly 
larger C0 2 amplitude than seen in Figure 25, and at least part of the deficiency around 
the tropopause is due to the coarse vertical resolution (2 km) of the model. Because 
of this resolution, we must ramp down I\ zz to very small values starting in the very 
upper troposphere (~15 km) to avoid spurrious transport across the tropical tropopause 
(section A.2.b). In reality, this region experiences large vertical mixing, and so the 
seasonal C0 2 signal is propagated through the tropopause with a larger amplitude than 
is simulated by the model. Future versions of the model with improved resolution should 
provide better simulations of the region around the tropopause. 

To check the propagation of the C0 2 seasonal cycle throughout the lower 
stratosphere, Figure 26 compares the model simulation (bottom row) with time series 
representing the ER-2 data at several latitudes and levels adapted from Strahan et al. 
[1998] (top row). These curves are least squares fits to a combination of the estimated 
stratospheric input seasonal cycle (the delayed Mauna Loa-Samoa average in Figure 25) 
and the observed linear trend at stratospheric midlatitudes. The data are binned by 
N 2 0 for the range of potential temperature (0) values indicated (370-410K, 395-435K, 
440-480K), and for three latitude ranges: tropics (6°S-12°N), subtropics (12°S-30°N), 
and midlatitudes (30°S-48°N). The model values represent the averages for the same 
latitude and theta ranges used for the observations. In all three latitude regions, the 
observations show similar amplitudes in the lower two theta ranges. At the top level, the 
data show some amplitude attenuation in the tropics, with much more attenuation with 
increasing latitude so that little seasonal signal is seen at midlatitudes at 440-480K. The 
data also show a general increasing phase lag with increasing altitude and latitude. 

The model simulation in Figure 26 shows an underestimation of amplitude at all 
levels and latitudes. This is due in large part to the underestimation of the signal 
amplitude entering the stratosphere (Figure 25) which then biases the model simulation 
throughout the stratosphere. Nevertheless, the simulation shows some qualitative 
consistency with the observations and illustrates that the model reproduces the general 
characteristics of seasonal cycle propagation. For example in the tropics and subtropics, 
the model and data show an increasing phase lag and amplitude attenuation with 
height. At midlatitudes, the model shows very good qualitative agreement with the 
observations: very similar amplitudes and phases are seen in the bottom two levels, with 
the top level showing almost no seasonal amplitude. At the lowest level (370-410K), the 
model compares well with the observations in simulating very similar amplitudes and 
phases in the tropics and subtropics (Figure 27), illustrating strong coupling of these 
regions. The midlatitude data at 370-410K show some amplitude attenuation and phase 
delay relative to the tropics and subtropics, and these features are reasonably captured 
by the model, although more amplitude attenuation is simulated than observed. 

Referring back to Figure 23a, somewhat surprising features are the simulated 
amplitude maxima at middle and high latitudes of both hemispheres centered near 16 
km. At these altitudes and latitudes, the seasonal cycle due to the C0 2 source has been 
completely attenuated. However, descent within the polar vorticies during winter and 
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early spring advects down relatively low C0 2 values, creating a late winter /early spring 
minimum and hence a significant seasonal cycle amplitude. This effect maximizes at the 
poles, and is evident at 15-20 km where there is a substantial vertical gradient in C0 2 . 
Above 20 km, there is little vertical gradient in C0 2 so that the downward motion has 
little effect. This feature was not indicated in previous model simulations of the C0 2 
seasonal cycle [Hall and Prather, 1993; Waugh et al., 1997; Strahan et al., 1998]. This 
feature would also be difficult to isolate in the ER-2 data given the limited temporal 
coverage of the current observations. Furthermore, this feature does not seem to be 
related to the annual cycle in the C0 2 source, so that extracting information from the 
data by fitting a tropospheric annual cycle and trend as was done by Strahan et al. 
[1998] would not seem to be justified in this case. 

Figures 28 and 29 compare vertical profiles of the ER-2 and model C0 2 for different 
months in the tropics, northern midlatitudes, and high southern latitudes. Here, C0 2 
is plotted against potential temperature, which acts as a vertical coordinate similar to 
pressure. The secular trend in C0 2 is apparent in the observations, and is captured fairly 
well by the model (e.g., Figures 28a,f). The tropical data in Figure 28 show the vertical 
propagation of the C0 2 seasonal cycle [Boering et al., 1996] which is repeated from 
year to year. A maximum is observed just above the tropical tropopause (380-400K) 
during August, and progresses to higher levels during October through December. A 
minimum is seen near 380K during October and at 430-440K during February. The 
model shows some indication of this signal propagation, however, the amplitude is 
significantly underestimated, as was seen in Figure 26. As shown in Figure 25, the 
model underestimates the C0 2 seasonal signal entering the stratosphere by roughly a 
third. However, Figure 28 suggests that relative to the observations, the model further 
overattenuates the amplitude with increasing height in the tropical lower stratosphere. 
A similar deficiency was seen in the propagation of the seasonal cycle of 2CH 4 + H 2 0 
(Figure 9). 

At middle and high latitudes (Figure 29), the vertical variations of the model show 
general agreement with the observations. Both the model and data indicate a larger 
decrease with height at midlatitudes compared to either the tropics or high latitudes. 
The large amount of scatter seen in the observations of C0 2 and 0 in Figure 29 can 
be explained by the fact that in the extratropics, air parcels originate from a wide 
range of airmass types due to the inherent large dynamical variability. And unlike C0 2 , 
potential temperature is conserved for time scales of only 2-3 weeks at most in the lower 
stratosphere. Therefore, C0 2 and 0 do not exhibit a high correlation in these regions. 
Furthermore, these data reflect variations in transport for a particular year that will 
not be resolved by our climatological model. However in the tropics, air parcels on 
a given theta (or pressure) surface have a limited distribution of air mass origins, so 
that C0 2 and 0 form a compact relationship [e.g., Boering et al., 1996]. To reduce the 
variability due to transport in the extratropics, it is useful to plot C0 2 using N 2 0 as the 
vertical coordinate [Boering et ah, 1996; Strahan et al, 1998], as shown in Figure 30. 

In the lower stratosphere, both C0 2 and N 2 0 are very long lived and form a compact 
relationship, unlike C0 2 and 0. The tropospheric source of N 2 0 is seasonally invariant 
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and undergoes only a small secular increase which is contained in our model simulation. 
Therefore, variations in CO 2 for a given value of N 2 0 will be independent of the seasonal 
variations in transport, and will be due only to the seasonal cycle and trend of the C0 2 
source. 

The model C0 2 — N 2 0 relationship shows good overall consistency with the data. 
Note the observed C0 2 maximum at N 2 0 = 300 ppbv during October/November 1995 
is qualitatively simulated by the model, although with reduced magnitude. This feature 
probably reflects the separation between tropospheric air at the lowest altitudes (highest 
N 2 0 concentrations), and stratospheric air above. A notable discrepancy in Figure 30 
appears at the high southern latitudes where the model consistently underestimates the 
vertical C0 2 gradient at low N 2 0 values. This results in too much C0 2 being simulated 
by the model at these upper levels, a deficiency which occurs throughout the year. At 
these altitudes and latitudes, the seasonal cycle is almost completely attenuated so that 
C0 2 exhibits only a secular increase. The model-measurement difference may therefore 
be due to inadequacies in the model photochemistry such that, 1) the N 2 0 loss rates 
are too large; 2) the production of C0 2 from CH 4 oxidation is overestimated, and/or 3) 
there is an additional stratospheric loss of C0 2 not accounted for in the model. 

4* Conclusions 

We have developed a new empirically-based transport algorithm for use in our 
GSFC two-dimensional transport and chemistry model. To derive the transport fields, 
the algorithm utilizes empirically based zonal mean zonal wind, temperature, diabatic 
and latent heating rates, E-P flux diagnostics for planetary and synoptic scale waves, 
and parameterizations to account for the effects due to gravity waves and equatorial 
Kelvin waves. As such, this scheme utilizes significantly more information compared to 
our previous algorithm which was based only on zonal mean temperatures and heating 
rates. 

To validate this new transport, we have presented extensive model-measurement 
comparisons of several long lived tracers. The model captures much of the qualitative 
spatial and seasonal variability seen in satellite and aircraft observations. These include 
the isolation of the tropics and the well mixed surf-zone region of the winter sub-tropics 
and mid-latitudes. The fact that the simulated tracer gradients in the tropics and 
subtropics agree fairly well with the observations supports the validity of using the 
HRDI wind data to derive the model wave driving and diffusion quantities in this region. 
At high latitudes, the model captures the general isolation of the southern polar vortex 
in winter and early spring, but does not adequately resolve the NH polar vortex. In the 
tropical upper stratosphere, the simulated double peaked structure in CH4 associated 
with the SAO is similar in phase, but is significantly weaker in amplitude compared to 
observations. This deficiency is consistent with other previous model studies. 

The model total ozone field captures much of the observed seasonal and latitudinal 
variability seen in the TOMS data. A major discrepancy occurs in the middle to 
high southern latitudes in which the model does not resolve the spring break-up of 
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the Antarctic ozone hole and subsequently maintains larger than observed latitudinal 
gradients during the southern high latitude summer. Sensitivity tests reveal that this 
deficiency is not due to an underestimation of horizontal mixing during the break-up of 
the southern polar vortex, and points to other processes that are not accurately resolved 
in the model. 

In the mesosphere, the model simulated H2O also shows good agreement with 
HALOE data. The model reveals that the seasonal water vapor changes are controlled 
primarily by vertical advection by the meridional circulation. Vertical diffusion induced 
by gravity wave breaking is less important, consistent with the findings of Holton and 
Schoeberl [1988]. However, diffusion does generate a secondary H 2 0 maximum during 
mid-winter, which is somewhat more pronounced in the southern hemisphere. 

We performed several diagnostic tests of the model transport emphasizing the 
region below 35 km. Simulations of carbon- 14 and strontium-90 reproduced the 
observed peak in mixing ratio at 20-25 km, and the decrease with altitude in mixing 
ratio above 25 km. The simulated mean age from SFg reveals that the oldest air of 
around 6 years occurs in the high latitude upper stratosphere during fall and early 
winter. The latitudinal and vertical gradients of the simulated mean ages compare well 
with aircraft and balloon observations outside of the northern hemisphere polar winter 
vortex. Including a mesospheric loss rate identical to that of CFC-115 (corresponding 
to a ~600 year atmospheric lifetime) in the SF6 simulation increased the age of air by 
as much as 3 years in the high latitude middle stratosphere, and provided better overall 
agreement with the balloon data above 20 km. 

The fact that the model significantly underestimates (by 3-5 years above 25 km) the 
NH vortex observations of very old air of 8-10 years seems to be inconsistent with the 
reasonably good model-measurement agreement at other latitudes. Even when including 
the mesospheric loss for SFg, the model mean age at high northern winter latitudes still 
underestimates the observations by as much as 2 years. Some underestimation of the 
model 7 in this region is expected given the model underestimation of the isolation of 
the NH vortex seen in the CH4 comparison (Figure 2). The model does qualitatively 
show the greater isolation of the SH vortex, as indicated in the various simulations and 
the fact that the mean age is .75 years older in the SH than the NH at 20 km. If the 
observations indicating mean ages of 8-10 years within the NH vortex are accurate, the 
significant model underestimation of 7 in this region probably reflects the limitations 
of the current zonal mean climatological model formulation. Such a model would not 
be able to resolve the large longitudinal and interannual variability characteristic of the 
NH polar regions during winter and early spring. 

The propagation of the seasonal signal in carbon dioxide from the tropical 
lower stratosphere to midlatitudes is qualitatively reproduced by the model, and the 
model simulates the observed near-absense of seasonal variability at higher levels at 
midlatitudes. However, simulations of both C0 2 and H show that the model attenuates 
the seasonal signals too rapidly with altitude in the tropics, a characteristic that is 
consistent among many stratospheric models. We found that the model reproduces 
the observed altitudinal variation of horizontal entrainment of midlatitude air into 
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the tropics reported in several studies [e.g., Schoeberl et al , 1997; Mote et al f 1998]. 
Furthermore, the vertical phase propagation of H and the vertical gradient of the 
mean age are both consistent with the observations. These results, along with various 
sensitivity tests, suggest that the overattenuation of the seasonal signal cannot be 
explained by an over-estimation of the horizontal in-mixing, or the explicit vertical 
diffusion contained in the model. The numerical advection scheme exhibits small 
numerical diffusion [Lin and Rood, 1996], so that this would have only a small effect on 
the signal degradation. It is possible that the signal over-attenuation is caused in large 
part by the coarse vertical and horizontal model resolution (2 km by 10 ° latitude). The 
model with such resolution cannot sufficiently resolve and propagate the seasonal cycles 
so that the amplitudes are effectively reduced. 

We found that the circulation in the upper troposphere and stratosphere is 
determined by a combination of the empirically determined heating rates, wave 
drive, and the bottom boundary condition in the lower troposphere. Changes in this 
circulation will affect methane, total ozone, and the mean age in a consistent manner. 
For example, an overly strong circulation will overestimate (underestimate) CH4 at low 
(high) latitudes, with the opposite tendency in total ozone, and produce a mean age 
that is too young at all latitudes. Obtaining reasonable model-measurement agreement 
in these tracers simultaneously at different latitudes and seasons gives confidence 
that the circulation is fairly consistent with the real atmosphere. Furthermore, the 
simulations with the new transport formulation are in substantially better agreement 
with observations compared to our previous model transport. Taken together, these 
model-measurement tracer comparisons therefore provide a rigorous test of the new 
model transport, and enable us to provide better model ozone simulations for a variety 
of scientific and assessment studies. 
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Appendix: GSFC 2-D Model Transport Description 

Earlier versions of the GSFC 2-D model transport algorithm derived the residual 
meridional and vertical velocities from the TEM thermodynamic and continuity 

equations. This formulation used empirically determined climatological temperatures and 
heating rates and neglected the eddy heating terms. This was essentially the diabatic 
circulation, as based on Dunkerton [1978], and assumed that horizontal temperature advection 
by the circulation was small. This assumption is reasonable in the stratosphere, but may not 
be applicable in the mesosphere where meridional velocities can be large. This approach was 
somewhat different from other studies that derived the diabatic circulation by iterating the 
continuity and thermodynamic equations to obtain a converged solution [Shine, 1989; Huang 
and Smith, 1991, Eluszkiewicz et ah, 1996]. In the previous GSFC model versions, the mean 
residual velocities along with the empirical zonal mean wind field were used in the zonal mean 
momomentum equation to obtain the Eliassen-Palm (E-P) flux divergence, or equivalently, 
the eddy flux of potential vorticity (PV) taken as the residual needed for momentum balance. 
The K yy field was then determined self-consistently as the ratio of the eddy PV flux to the 
latitudinal gradient of zonal mean PV [Newman et al., 1988; Fleming et al ., 1995]. These 
circulation and diffusion fields were able to capture some of the general features seen in 
long-lived tracer observations. However, because the transport was based only on zonal mean 
temperatures and heating rates with no explicit planetary or gravity wave effects included, 
many detailed observational features reflecting zonal mean transport were not resolved in the 
model simulated tracer fields. 

Our new methodology generally follows that of other interactive 2-D models [e.g., 
Brasseur et al., 1990; Garcia et al, 1992; Bacmeister et al., 1995] as originally formulated 
by Garcia and Solomon [1983], and includes explicit calculations of the various wave effects. 
However, the new dynamical fields are based on pre-determined empirical data sets, as opposed 
to being computed interactively in the model. Note that as with previous model versions, the 
transport is climatological and does not account for the quasi-biennial oscillation or any other 
interannual variability. 


A.l. Zonal Mean Formulation 


Following the methodology originally outlined in Garcia and Solomon [1983], a meridional 
streamfunction (x*) is calculated to obtain the TEM residual circulation (v*,u>*). To derive 
an equation for the streamfunction, we use the zonal mean momentum, thermodynamic, and 
mass continuity equations of the TEM system [e.g., Andrews et al., 1987], 


du utan<^ 

dt a ad<j> 


d u . du _ 

) + w — = Fpw + Fqw + Fkw 


dt m d T m/ dT kT. . 

W + ” rti + * ( Tz + = Qd + Ql + Qpw + Qaw 


1 dp* cos (ft) 1 djvfp) _ 
a cos <f> d(f> pdz 


(1) 

( 2 ) 
( 3 ) 


where the scale height ( H ) is taken to be 7 km, and the other symbols have their usual 
meanings. Here, mechanical forcing is taken as the E-P flux divergence due to dissipating 
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planetary waves (Fpw), gravity waves {Fqw)^ and equatorial Kelvin waves ( Fkw )* We have 
not included a term accounting for unresolved motions, e.g., as represented by a weak Rayleigh 
frictional drag, given the uncertainties involved. Thermal forcing is due to the combined effects 
of the net diabatic heating (Qd)i latent heating ( Qp ), and the net eddy heat flux associated 
with planetary waves (Qpw) and gravity wave breaking in the mesosphere ( Qgw )■ Equations 
(1) and (2) are combined using the thermal wind relation, 

du R dT 2u tan<£ 

/ 1 r\ — tt jO /l — / ' ' V^J 

oz H ao<p a 


This eliminates the time derivatives (assuming u z df\/dt is negligible), and along with 
continuity yields a single diagnostic equation for x*, 

d 2 \ * dx* d 2 x* dx* d 2 x* 

c --& +c >ik + + c + c «w + Cor = CFms * (5) 

The coefficients {C ZZ ,C Z , C zy , C y , C yy , and Co) depend on empirically determined u and 
T fields (see section A.l.b) and are similar to those derived in Garcia and Solomon [1983]. The 
forcing term ( Cp ) is proportional to the vertical gradient of the total momentum forcing and 
the latitudinal gradient of the total heating rate, 

d R d 

Cf = J\-Q^{Fpw + F C w + Fp-w) + + Ql + Qgw + Qpw ) (6) 


These quantities will be described below. The residual circulation is computed from x * via 
continuity, 


COS 4> 


1 


w = 


dx* 

8z " 

dr 


r' 

H 


cos <p dy 


(7) 

( 8 ) 


Note that this formulation satisfies global mass balance such that, 



pw * cos 4>d<f) = 0 


(9) 


This differs from previous calculations of the diabatic circulation in which inaccuracies in the 
heating rate algorithms and/or input data sets made it necessary to manually adjust w* to 
satisfy the condition in (9). 


A. l.a. Numerical Solution/Boundary Conditions on x * 
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To obtain x*> equation (5) was solved in finite difference form by a successive over- 
relaxation method (identical results were obtained using the method of Lindzen and Kuo 
[1969]). We used a rectangular latitude-height grid extending from 90°S-90°N at 5° increments, 
and from the ground to ~116 km in log-pressure altitude at approximately 2 km intervals. As 
discussed below, the top level of the transport model has been raised from the 90 km level 
used in previous model versions. 

For boundary conditions, we treat the model as being bounded by four rigid walls with no 
flow into or out of the model domain. Therefore, at the poles and u>*=0 at the ground 

(1013 mb) and the top boundary. By equations (7) and (8), this implies that x*=0 at the 
four boundaries. While our current chemistry model extends up to 90 km, putting the upper 
boundary on the dynamics and hence making w*=0 at this level caused unrealistically large 
values in v* in the upper mesosphere. To avoid this problem, we extended the upper boundary 
on the model transport up to 116 km. The heating rates and wave driving were linearly 
damped down to be zero at this top level. This gave a w* field that gradually decreased to 
zero at the top boundary, thereby avoiding any anomolously large jets in the v* field. In this 
formulation, the dynamical processes in the region above 90 km are not completely represented 
and as such this region serves only as a sponge layer for the transport. We therefore only show 
model results for the 0-90 km region where chemistry is computed. 

To obtain reasonable long lived tracer simulations in the stratosphere, we found it 
necessary to impose an additional boundary condition on x* at one level above the ground 
(762 mb - denoted as X2*)- Following Garcia et al [1992], we calculate X2* as, 

cos d) 

X2* = — 7 — / p{Fpw + Fgw + Fkw)<Iz (10) 

J P2 J 


This expression is derived from equation (7) and the linearized steady state momentum 
equation, — fv* = Fpw + Fqw + Fh 'W? and follows from the downward control principle under 
quasi-geostrophic scaling [e.g., Haynes et al. } 1991; Rosenlof and Holton, 1993]. Equation (10) 
is applied to latitudes outside of ±15° with values in the tropics interpolated between ±15°. 
The vertical velocity at 762 mbar obtained from X2* then depends solely on the momentum 
forcing of the free atmosphere above this pressure surface. Upward motion out of the boundary 
layer generally occurs equatorward of ±50° with downward motion into the boundary layer at 
higher latitudes. We found that the circulation computed using \2* obtained from (10) was 
too strong, i.e., too much upwelling (downwelling) in the tropics (high latitudes), as reflected 
in model-measurement comparisons of CH4, total ozone, and the mean age of air (too young at 
most latitudes). Multiplying the expression in (10) by a constant factor of .75 for all latitudes 
and seasons produced a somewhat weaker circulation and gave much better model agreement 
with long lived tracer observations. We therefore used this methodology for all simulations 
shown throughout this paper. We note that as an independent check of this method, a very 
similar circulation was obtained by computing X2* via integration of equation (8) and applying 
the definition [e.g., Andrews et al. , 1987], 
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Here, the terms on the RHS of (11) are obtained from climatological averages of the recently 
reanalyzed NCEP data [Kalnay et aL, 1996]. 

Figure A.lc-d shows w* and v* for January. These fields are generally similar to previous 
model calculations [e.g., Garcia et al. } 1992]. Upward motion is seen throughout the tropical 
stratosphere, and in the extratropical middle to upper stratosphere and mesosphere in summer. 
Maximum ascent of 3.6 cms~ l occurs in the polar upper mesosphere. The winter extratropics 
are characterized by downward motion throughout the troposphere and middle atmosphere. 
The strong upward motion of the Hadley cell is also evident in the tropical troposphere. The 
meridional velocities are generally weak in the stratosphere, but reach maximum values of 
6 — 7 ras -1 in the upper mesosphere reflecting the summer to winter meridional flow pattern 
induced by gravity wave drag. 


A. Lb. u and T 

The zonal mean temperature and zonal wind fields used in the coefficients in equation (5) 
were compiled as follows. We used zonal mean temperatures based on the 17-year (1979-1995) 
average of NCEP data for 1000-1 mbar, and the CIRA-86 empirical reference atmosphere 
[Fleming et aL , 1990] for the mesosphere above 1 mbar. Outside the tropics, zonal mean winds 
were derived from the temperatures via thermal wind balance (equation 4). This derivation is 
problematic at low latitudes where the Coriolis parameter is small, and in the region above 80 
km where the temperature data contained in the CIRA model was rather sparse. Therefore, 
we used direct wind measurements from the high resolution Doppler imager (HRDI) onboard 
UARS to obtain u for 20°S-20°N in the stratosphere and mesosphere (10-115 km) and at all 
latitudes above 80 km. We used the HRDI level 2B winds averaged over the period November 
1991 to September 1996 which cover the region 70°S-70°N. To maintain a constant relative 
angular velocity at high latitudes, values of u poleward of ±70° were extrapolated by cos <j> so 
that u — > 0 at the poles. The HRDI data also contain a measurement gap at 40-50 km which 
was filled in by spline interpolation. Above ~75 km, the HRDI data at lower latitudes (±35°) 
becomes increasingly aliased by the diurnal tide due to the daytime only measurements. This 
aliasing was minimized by fitting and subtracting out the tidal mode corresponding to the (1,1) 
Hough function, expected to be the dominant contribution. This fitting procedure is described 
by Hays et al. [1994] and Burrage et aL [1995]. Since HRDI does not make measurements 
below 10 km, we derived the zonal wind in the equatorial lower-middle troposphere from 
a second derivative thermal wind calculation. In this formulation, u is proportional to 
( df/dy ) _1 , rather than / _1 , thereby avoiding the equatorial singularity problematic of the 
standard geostrophic wind derivation. This technique has been used in the middle atmosphere 
with some qualitative success in estimating the directly measured wind [e.g., Fleming and 
Chandra, 1989; Fleming et al. , 1996]. 

As noted above, the streamfunction formulation necessitates maintaining exact thermal 
wind balance in all regions. To do this in regions where different temperature and wind data 
were used, we combined the du/dz field computed directly from the HRDI winds with the 
du/dz field implied from the temperatures via thermal wind balance. These quantities were 
merged together over several latitudes or levels to obtain a smooth transition. Again, we set 
the constraint that at the poles, u = 0 and therefore, du/dz — dT/dy — 0. The resulting 
du/dz and the equivalent dT/dy fields were then re-integrated to obtain u and T fields in 
thermal wind balance. Figure A.la-b shows these fields for January and indicates the well 
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known zonal mean structure during solstice conditions, e.g., the summer easterlies and winter 
westerlies of the middle atmosphere. These features have been discussed frequently in previous 
climatological compilations [e.g., Randel, 1992; Fleming et al, 1990]. 


A.l.c. Heating Rates 

The diabatic heating rates (Qd) are computed following Rosenfield et al. [1994], using 
climatological distributions of zonal mean temperature (section A.l.b), ozone, and water 
vapor. The ozone climatology was constructed as follows. Values in the troposphere are based 
on SBUV data as described in McPeters et al [1984], SAGE II data (both sunrise and sunset) 
averaged over 1984-1989 are used in the lower stratosphere (100-10 mbar). Nimbus-7 SBUV 
data averaged over 1979-1989 are used for the upper stratosphere (10-0.7 mbar), and SME 
IR and UV measurements averaged over 1982-1986 are used for the mesosphere above 0.7 
mbar. The water vapor climatology is based on Oort [1983] for the troposphere up to 375 
mbar and climatologically averaged SAGE II data for the upper troposphere and stratosphere 
(375-1 mbar). For the mesosphere above 1 mbar, we used the midlatitude reference profile 
from Remsberg et al [1989] which is independent of latitude and season. For both the ozone 
and water vapor climatologies, linear interpolation was used to fill in data gaps. For data-void 
regions at polar latitudes, data were extrapolated by maintaining a constant value poleward 
from the highest latitude available. Values were blended over 2-3 levels where there was 
overlapping data, thereby obtaining a smooth transition between the original data sources. 

The net diabatic heating rates for January are shown in Figure A. 2. These values 
are generally consistent with previous heating rate calculations that utilize climatological 
data [e.g., Shine, 1989; Huang and Smith, 1991; Newman and Rosenfield, 1997], and UARS 
observations [Eluszkiewicz et al, 1996]. Net diabatic heating (cooling) is seen in regions where 
the observed temperature is colder (warmer) than the radiative equilibrium temperature. 
Largest cooling occurs in the winter extratropics due to planetary and gravity wave effects, 
with the latter being dominant in the mesosphere. Gravity wave breaking drives the summer 
mesosphere to be far colder than the radiative equilibrium temperature, leading to large net 
diabatic heating there. Small net diabatic cooling occurs throughout the troposphere during 
all seasons. Largest values are ~-1.5 K/day which occur in the tropical mid- troposphere. 

For the latent heating rates ( Ql ), we have used the empirical climatology compiled by 
Newell et al [1974]. The January latent heating is shown in Figure A. 3a (note that this plot 
is only for the troposphere, 1000-100 mb). Largest latent heating of ~2.5 K/day occurs in the 
tropical mid-troposphere coincident with tropical convective systems, with minima observed 
in the high latitude regions. 

Gravity wave breaking in the mesosphere induces a net thermal flux due to the 
combination of an eddy heat flux convergence and the turbulent diffusion of the mean thermal 
field [Schoeberl et al, 1983; Huang and Smith, 1991]. The time tendency due to the eddy heat 
flux convergence is expressed as, 
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where k = R/C p , c is the gravity wave phase speed, a is the Newtonian cooling rate 
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associated with the radiative damping of the waves, and K zz is the vertical eddy diffusion 
coefficient computed from the gravity wave breaking paramerization described in section 
A.2.b. Consistent with Huang and Smith [1991], our calculations show that the first term in 
parentheses in equation (12) is quite small compared to the turbulent diffusion term. The total 
quantity, Qgw i in (12), is shown in Figure A. 3b for January. Net cooling is seen throughout 
the mesosphere, with largest values of ~-4.5 K/day occurring in the regions of largest diffusion 
in the upper mesosphere (see Figure A.5d). 

Turbulent diffusion of the mean thermal field is expressed as, 
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This quantity is shown in Figure A. 3c for January and can contribute to either a net heating 
or cooling depending on the vertical distribution of temperature and diffusivity. For January, 
net cooling is seen throughout most of the mesosphere with a maximum of -11 K/day in 
the summer mid-latitudes near 80 km. Smaller regions of net heating on the order of +1-2 
K/day occur in the middle mesosphere at mid-latitudes of both hemispheres. The total 
contribution of breaking gravity waves to the thermodynamic budget (equation 2) and used 
in the streamfunction (equations 5-6) is then Qgw = Qgw\ + QgW 2 shown in Figure A. 3d. 
The total net contribution of gravity wave breaking is cooling throughout almost the entire 
mesosphere. A maximum of -12.5 K/day occurs at 80-90 km in the mid-latitude SH during 
January. 

Planetary waves also induce a net thermal flux [e.g., Andrews et al. } 1987] which enters 
into the streamfunction equations (5)-(6) as, 
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The eddy terms in (14) are taken from 3-D meteorological fields as discussed in section 
A. 2. a. The quantity Qpw (figure not shown) is small compared to the other heating 
rate contributions. Maximum values are generally less than ±.75 K/day occurring in the 
high latitude winter stratosphere-lower mesosphere, and the mid-high latitude troposphere 
throughout the year. 

The total heating rate for January is shown in Figure A. 4. Diabatic heating/cooling 
dominates the total throughout the stratosphere, with the gravity wave effects playing a 
significant role in the mesosphere along with the diabatic contribution. Both the diabatic and 
latent heating are important in the troposphere, and the combination of these processes leads 
to horizontal gradients in the total heating rate which determines the classical Hadley cell 
circulation in the tropical troposphere and lower stratosphere (Figure A.lc). 

A. 2. Wave Parameterizations 

A. 2m. Planetary and Synoptic Scale Waves 

Planetary wave dissipation generates rapid irreversible mixing of tracers which can be 
expressed as a horizontal diffusive transport in two-dimensional models. This wave dissipation 
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also induces a drag on the zonal mean flow which can modify the mean advective transport 
by the meridional circulation. Previous investigations have illustrated the importance of 
accounting for these processes in a self-consistent manner [Newman et ai, 1988; Garcia, 1991], 
The E-P flux divergence computed from 3-dimensional meteorological analyses provides a 
diagnostic estimate of the planetary wave drag ( Fp\y in equation 1). This quantity is defined 
as [e.g., Andrews et al. , 1987], 
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Following Randel and Garcia [1994], horizontal mixing rates can then be calculated self- 
consistently as the ratio of the E-P flux divergence to the latitudinal gradient of zonal mean 
potential vorticity, 

Kyy — - Fp\V/Qy ( 18 ' 
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This formulation follows from the flux-gradient relationship [e.g., Newman et al. , 1988], 
v f q f = — Kyy q y , where the E-P flux divergence is taken as the zonal mean PV flux for 
quasi-geostrophic planetary waves. 

For the troposphere and stratosphere, we have derived climatological E-P flux divergence 
and Kyy quantities from equations (15)-(19) using 17 years (1979-1995) of daily 3-dimensional 
NCEP analyses for 1000-1 mb. This computation includes the wave dissipation and induced 
mixing due to both thermal damping and wave breaking. To account for the effects from both 
planetary and synoptic scale waves in the troposphere and very lower stratosphere, we include 
zonal wavenumbers 1-12 from the ground to 10 mb. Above this level, zonal waves 1-6 are 
included since only the planetary scale waves are important in this region. 

For the zonal and meridional velocities in ( 15)-( 19), balanced winds were computed from 
the NCEP geopotential height analyses for regions poleward of ±20° [Randel, 1992]. To deal 
with the problem of deriving winds in the tropics, we computed balanced winds at the equator 
from the NCEP height fields using a second derivative calculation analogous to that used for 
the zonal mean wind as described in section A. l.b. Winds between ±20° and the equator were 
interpolated by the cosine of latitude to minimize the latitudinal gradient across the equator 
and thereby maintaining inertial stability. Previous studies have shown that this techique 
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produces a u field that is in general qualitative agreement with radiosonde, rocketsonde and 
HRDI data in the phasing of the QBO and SAO. However, the derived wind magnitudes 
were generally underestimated compared with the direct wind measurements [Fleming and 
Chandra, 1989; Fleming et al, 1996]. We found that deriving the equatorial eddy winds in 
this manner underestimated the E-P flux and eddy diffusion magnitudes. This was evident in 
the simulated long lived tracer fields that exhibited significantly larger latitudinal gradients in 
the tropics compared to UARS observations. We also found that using tropical eddy winds 
from the UKMO data assimilation system provided only a marginal improvement in the 
simulated tracer fields. Instead, we utilized HRDI zonal and meridional wind measurements 
for the calculations between ±20°. To avoid a bias caused by the QBO, we used HRDI data 
for January 1992 to May 1996 corresponding to 2 complete QBO cycles. Because of data 
gaps and the limited longitudinal coverage obtained with daily HRDI data, we used monthly 
mean HRDI winds and the corresponding monthly mean NCEP temperatures to compute 
the quantities in equations (15)-(19) between ±20°. Using monthly mean data accounts for 
the stationary wave components but not the transient eddy contributions. However, the 
HRDTderived E-P flux divergence and K yy values were larger than those derived from the 
NCEP or UKMO winds. The resulting tropical-subtropical gradients in the tracer fields were 
weaker and in significantly better agreement with observations compared with the NCEP or 
UKMO-based derivations (see Figures 1 and 19). 

For the mesosphere above 1 mbar (up to ~ 85km), we compute the E-P flux divergence 
and eddy mixing from the CIRA-86 stationary planetary wave climatology [Barnett and 
Labitzke, 1990]. This contained monthly mean geopotential height and temperature values for 
zonal waves 1-2 which are based on data from the Nimbus-6 Pressure Modulator Radiometer 
(PMR) for 1975-1978. As in the stratosphere, balanced winds were used for the wind velocities 
poleward of ±20°, with the second derivative balanced wind calculation used in the tropics. 
Unlike the stratosphere, we did not use the HRDI data for the tropical winds in the mesosphere 
since using different time periods for the temperatures (1975-1978) and winds (1992-1996) 
was thought to be unsatisfactory for the calculations in equations (15)-(19). Values of Fpw 
above 85 km were linearly damped down with altitude to be zero at 116 km. Future E-P flux 
calculations for the mesosphere will use both wind and temperature data from UARS as these 
retrieval algorithms mature. 

Consistent with previous investigations [Robinson, 1986; Newman et al. , 1988], our 
computations of equations (15)-(17) revealed mostly negative areas of E-P flux divergence, but 
also some positive values, especially at high latitudes. This implies a wave-induced acceleration 
of the mean flow, and a negative mixing rate via the flux-gradient relationship in (18). While 
regions of positive values may be caused by limitations of the original data sets, some areas 
of positive E-P flux may be real, and the implied negative K yy indicates a breakdown in the 
assumptions underlying (18). This is also indicative of the ambiguity in interpreting E-P 
diagnostics for wave propagation and dissipation [e.g., Andrews, 1987; Newman et al, 1988]. 

To avoid computing negative mixing rates in practice, regions of positive E-P flux 
divergence were set to zero, implying a zero K yy . The resultant E-P flux fields were smoothed 
once in the latitude-height domain before being applied in the calculations for x* and K yy . 

To ensure that K yy computed from (18) remains well behaved, we also apply the criterion 
that q y > 0.5xl0 _11 m _1 s _1 as suggested by Randel and Garcia [1994]. This avoids negative 
mixing rates and implies that large positive values of K yy are caused by eddy processes via 
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the E-P flux divergence and not by instability of the zonal mean flow as represented by small 
or negative q y . To accomodate the one day time step used in the chemistry model, we set 
an upper K yy limit of 5xl0 :o cm 2 s~ l which is sufficiently large to accomodate the diffusion 
induced by planetary and synoptic scale waves. 

We impose a very small background K yy of 10 8 cm 2 ^ 1 throughout the stratosphere and 
mesosphere. In the low-middle troposphere, it was necessary to set a large background K yy of 
2xl0 lo cra 2 s” 1 to simulate the interhemispheric transit time of 1-2 years estimated from SF6 
and halocarbon data near the ground [e.g., Geller et al, 1997]. This background K yy was 
gradually tapered down with altitude to blend with the small 10 s cm 2 s~ 1 background value 
imposed in the stratosphere and mesosphere. Finally, we impose a small K yy of 10 8 cra 2 s _1 
around the tropopause region where the tropopause height changes across adjacent latitudes. 
This decreased the anomolous cross-tropopause transport and gave better model simulations 
of water vapor compared with observations in the very lower stratosphere at middle to high 
latitudes. 

We do not include the transport due to the off-diagonal K yz and K zy terms in the 
diffusion tensor which arise from the projection of K yy from isentropes to pressure surfaces 
[Newman et al, 1988]. This was done because of the uncertainties in the calculation of K yz 
and K zyy and because these terms represent only a small effect. Note also that to ensure that 
the diffusion matrix is always positive and diffusive, the relation K yy K zz > K yz K zy must hold. 
Neglecting the off-diagonal terms ensures that the total diffusive flux is positive even when 
imposing the very small K zz values necessary in the stratosphere. 

Figures A. 5a and A. 5c show latitude-height plots of the climatological E-P flux divergence 
( Fpw ) and K yy for January. In the troposphere, large wave driving is evident poleward of 
±20° during winter and to a lesser extent, summer, due to synoptic and planetary wave 
activity. These features are similar in both hemispheres (July not shown). The K yy values 
computed from equation (18) reflect this pattern, although they are superimposed on the 
large background mixing rates which we have specified throughout the troposphere as 
described above. Large planetary wave driving also occurs in the extra-tropical stratosphere 
and mesosphere during winter. Consistent with Randel and Garcia [1994], the direct ratio 
calculation used in (18) gives large mixing rates at high latitudes in the upper stratosphere and 
lower mesophere, and also above 100 mb at 20°-45°N coincident with weak q y in the surf-zone 
region. This is in contrast to Newman et al [1988] who computed largest mixing rates only at 
high latitudes using a linear regression fit of v f q f and q y with a zero intercept. 

The seasonal variation of K yy and u at several different levels are shown in Figure A.6a-d 
(see also Figs. 12-13 in Randel and Garcia [1994]). These plots depict the hemispheric 
asymmetries in mixing and the zonal mean flow. In the mid-upper stratosphere and 
mesosphere, weak mixing occurs in the easterly flow regimes of the summer hemisphere and in 
the tropics, consistent with a minimum of planetary wave activity in these regions. Note that 
the tropical stratospheric mixing rates based on the ERDI winds are weak but not negligible, 
with values as large as 8xl0 9 cm 2 s~ 1 . In the mid-late winter, large K yy values are seen along 
the subtropical edges of the polar night jets, characterizing the surf-zone region of the winter 
stratosphere. This area of maxima moves poleward during late winter, with large mixing at 
high latitudes during spring of both hemispheres associated with the deceleration and break 
up of the polar vorticies. These figures also show the contrasting degree of winter polar vortex 
isolation in the two hemispheres. Weak mixing occurs in the core of the SH mid-winter vortex, 
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consistent with the inhibition of planetary wave propagation through strong westerly winds 
[Charney and Drazin, 1961]. The NH exhibits significantly weaker u and larger K yy throughout 
the winter stratosphere poleward of 25°N. This hemispheric asymmetry is also reflected 
in the long lived constituent simulations discussed in section 3. In the lower mesosphere, 
both hemispheres show strong subtropical jets near ±35° in mid-winter, with strong mixing 
occurring mainly on the poleward flanks of the jets. In the very lower stratosphere (68 mb), 
westerly winds and significant mixing occur throughout much of the year in the extratropics 
of both hemispheres. There is also evidence of a semi-annual signal in the mixing in both 
hemispheres, with largest K yy in spring and fall coinciding with the build-up and decrease of 
the winter polar jets. As at higher levels, weaker mixing at 68 mbar is observed in the tropics, 
and in mid-winter in the cores of the polar jets, especially in the SH. 


A. 2.b. Gravity Waves 

Computations of vertical eddy diffusion (K zz ) and mechanical forcing (Fqw in equation 1) 
from gravity waves are based on the parameterization originally formulated by Lindzen [1981] 
and modified by Holton and Zhu [1984] (hereafter HZ84). In this scheme, waves propagate 
vertically and break when the amplitudes become large enough to cause convective overturning. 
The resulting turbulent diffusion then inhibits further amplitude growth with height, leading 
to a momentum flux convergence. This mechanical forcing can accelerate or decelerate the 
zonal mean flow depending on the phase speed of the wave relative to u. This in turn induces 
a meridional circulation, and, along with the turbulent diffusion, exerts a strong influence on 
the distribution of trace gases in the middle atmosphere. 

Other studies have used this parameterization to compute the time evolution of the 
zonal mean flow [Schoeberl et aL, 1983; HZ84; Garcia and Solomon, 1985]. However, since 
u and T are pre-determined in our model, we compute the seasonal and spatial distribution 
of wave drag and diffusion based on the empirically determined zonal wind and temperature 
fields and an input set of gravity wave parameters. These computed wave effects are highly 
dependent on the background wind field. For example, previous studies have observed greater 
wave drag during the SH winter compared with the NH [e.g., Shine, 1989]. It was postulated 
that the larger planetary wave amplitudes of the NH winter generate a greater range of wind 
speeds at different longitudes throughout the underlying stratosphere compared with the SH 
winter. This effectively filters out a greater range of the gravity wave spectrum, resulting 
in less gravity wave drag and diffusion in the NH winter mesosphere. To account for such 
longitudinal dependence, we use a 3-D (longitude-latitude-height) background zonal wind field 
for the calculations of wave drag and diffusion. The zonal mean of this u field is the same as 
that described in section A.l.b, with the longitudinal variability based on NCEP and HRDI 
data for 1992-1996 as discussed in section A. 2. a. The longitudinal grid resolution is 30° which, 
although rather coarse, should be sufficient to accomodate the planetary scale features of the 
middle atmosphere. Although the wind field is comprised of a 5-year average, the greater 
degree of zonal asymmetry during the NH winter compared to the SH is represented, and is 
reflected in the computed fields of wave drag and diffusion. 

In accordance with the tropospheric sources of gravity waves, we assume a spectrum of 
9 waves with zonal phase speeds c = 0, ±10, ±20, ±30, ±40 ms~ l . Following HZ84, we specify 
a horizontal wavelength of 200 km for these waves, consistent with the dominant horizontal 
wavelength of breaking waves observed in the mesosphere [Vincent and Reid, 1983]. Here we 
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have assumed that waves are only traveling zonally, although similar results were obtained 
by assuming an isotropic spectrum of waves as in HZ84. Garcia et al. [1992] obtained a 
reasonable simulation of the mesospheric zonal wind SAO in their model by including two 
additional tropical waves with phase speeds c = ±50 ms -1 and zonal wavelength of 800 km. 
We have followed this methodology and included these waves in our computations with the 
largest boundary forcing (see below) at the equator. 

For each latitude and longitude, the algorithm first requires determination of the breaking 
level Z}> r of each wave. Since more than one wave component in the spectrum can break at a 
given location, the z^ r of a given wave will be modified by the diffusion ( D ) induced by waves 
breaking at lower levels. The breaking level also depends on molecular diffusion ( D mo { ) in 
the lower thermosphere (as defined in HZ84), and on thermal damping of the wave which is 
parameterized by Newtonian cooling, so that 


%br = Z 0 + 
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where c and k are the wave phase speed and zonal wavenumber (corresponding to a 200 or 
800 km wavelength), a is the Newtonian cooling coefficient, and z 0 is the bottom boundary 
(375 mbar). Forcing at z 0 is determined by u — (woNo/k) 2 ^ 3 /\u 0 — c) 1 / 3 . Here, u?o is a vertical 
velocity that is specified as a function of the wave phase speed. As with the previous modeling 
studies, we assume that the small phase velocity waves have the largest bottom boundary 
forcing. Values of wq range from .03 ms~ l for the slow phase speed waves to .002 ms~ l for the 
fastest waves. We also skewed wq to have slightly larger values for the westerly waves, following 
Huang and Smith [1991], We impose only a weak geographical variation of u5o 5 with slightly 
larger values at midlatitudes compared with the tropics and polar regions. To mimic the 
generation of gravity waves by convection, we also specify slightly larger wq over subtropical 
land areas during summer following Wu and Waters [1996]. However, we note that in our 
current model formulation, the geographical and seasonal variability in the computed gravity 
wave effects is determined primarily by the background wind in the middle atmosphere. 

The turbulent diffusion D for each wave component is computed in the altitude range 
Zb r < z < z c , where z c is the critical level where c = n, such that, 
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Starting with the wave with the lowest z^ r (in which D = 0 in equation 20), the breaking levels 
and diffusion are recomputed for each subsequent wave, taking into account the total diffusion 
induced by waves breaking at lower levels. 

For altitudes in the range Zf, r < z < z c , the vertical momentum flux convergence for each 
component is, 
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Below z^ r , a particular wave does not generate diffusion, but undergoes damping from thermal 
dissipation, molecular diffusion, and the diffusion induced from other waves that are breaking. 
This wave damping results in a momentum flux convergence given by HZ84 (see their equation 
40) which we also include in our computations of F x . 

The computed wave drag and diffusion may exhibit discontinuities between the breaking 
levels of the different wave components. To smooth out these features and to allow for sporadic 
wave breaking below ^ r , we include an exponential decay of diffusion ( D sp ) and momentum 
deposition ( F xsp ) below z *> r for each wave. The diffusion and wave drag, including that due 
to sporadic wave breaking, are summed over all waves to obtain total distributions of these 
quantities (K zz and Fqw )• We also include the wave driving due to the turbulent diffusion of 
the zonal flow so that, 


Faw = <(F, + F,„) + ~ 


PK ; 


du 

~d~z_ 


and Ii zz = c(D + D sp ) (23) 


where e is an intermittancy factor discussed below. The K zz and Fqw values are then zonally 
averaged, and smoothed in latitude and height. 

As discussed in previous studies, the Lindzen parameterization gives unreasonably large 
values of diffusion and momentum deposition. This problem was alleviated by applying an 
intermitency factor (e) to the computed values in equations (21)-(23). This represents the fact 
that gravity wave breaking is not continuous, but rather is sporadic in time and space. For 
our calculations we found it necessary to set e = .025 for all wave components. The values of e 
along with wo used to determine z^ r were specified to compute a wave drag peaking at about 
100 ms~ l day~ l in the upper mesosphere as inferred previously [e.g., Holton, 1983; Fritts and 
Vincent, 1987]. This wave drag was in rough qualitative agreement but larger in magnitude 
compared with the mesospheric momentum residual obtained from the empirically determined 
diabatic circulation [e.g., Shine, 1989; Huang and Smith, 1991]. As discussed in section 3. 2. a, 
our computed values of wave drag and diffusion gave good agreement in the computed seasonal 
cycle of mesospheric H 2 O compared with UARS/HALOE measurements. It may be possible 
to further improve the simulations of mesospheric tracers by specifying the input parameters 
to have much stronger seasonal or spatial variability than is presently employed. However, we 
have not done this given the current lack of gravity wave observations, although recent work 
has begun to better characterize gravity waves on a global basis [e.g., Wu and Waters, 1996; 
Alexander and Rosenlof, 1996; Alexander, 1998]. 

We also include a contribution to the zonal mean Fqw and K zz from the breaking 
diurnal tide in the lower thermosphere as discussed by Lindzen [1981]. In accordance with 
UARS/HRDI observations [Lieberman and Hays, 1994] and theoretical estimates [Forbes, 
1982], we specify the largest values to be at the equator at 90-100 km, with a semiannual 
variation maximizing at the equinoxes. Values exponentially decrease in latitude and height. 
Largest values for March at 90 km at the equator are set to -10 ms" 1 day~ l for wave drag and 
30 m 2 s~ l for diffusion. 

The total zonal mean Fqw from gravity waves and the diurnal tide is then applied in 
equations (5)-(6). The total K zz including the background value (section A.2.d) is then used 
in the model transport of chemical constituents. We impose an upper limit of 100 m 2 s~ l to 
accomodate the 3 hour time step used for the vertical diffusive transport of constituents. 
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Previous studies have discussed the possibility that the turbulent Prandtl number, Pr 
(i.e., the ratio of the total eddy momentum diffusivity to the eddy diffusivity of heat or 
constituents) has a value greater than unity [e.g., Smith and Brasseur, 1991; Huang and Smith; 
1991]. Quantification of Pr in our model is difficult given that we do not compute u or T . 
Applying a value of Pr > 1 to the quantities computed in our gravity wave scheme essentially 
decreases the relative importance of the diffusive transport of constituents compared to the 
advective transport. However, the advective contribution already dominates in our current 
model formulation which has Pr = 1 [Chandra et ai, 1997]. Using a value of Pr > 1 did 
not improve the overall simulation of mesospheric H 2 O compared with UARS/HALOE data. 
Because of this and the uncertainties in determining Pr, we have elected to set Pr = 1 in the 
present model study. 

Figures A. 5b and A.5d show latitude-height sections of gravity wave drag (Few) and 
vertical diffusion from the gravity wave scheme together with the background specifications 
(see section A.2.d). These calculations are consistent with previous model results [e.g., Garcia 
and Solomon, 1985]. The largest wave drag occurs in the middle to upper mesosphere. Values 
of Fqw are opposite in sign to u (Figure A. lb) as the gravity wave drag acts to close off the 
jets in the upper mesosphere, and induce the summer-winter meridional flow pattern (Figure 
A. Id). The K zz field also shows maxima occurring in regions of large u in the mid-latitude 
upper mesosphere, with a distinct minimum in the tropics. This pattern is very similar to 
zonally averaged small scale variances seen in UARS microwave limb sounder (MLS) radiance 
measurements, which may be an indicator of gravity wave activity [Wu and Waters, 1996]. 
The tropopause is also evident in this figure, as seen by the sharp transition from large K zz in 
the troposphere to very small values in the stratosphere. The determination of the model Ii zz 
in this region is dicussed in section A.2.d. 


A.2.C. Kelvin Waves 

Equatorial Kelvin waves are thought to be important in driving the westerly acceleration 
of the stratospheric SAO and related meridional circulation. To approximate this process, 
we used the parameterization of Dunkerton [1979] and Gray and Pyle [1987] along with the 
empirical u and T (section A.l.b) to diagnose the mechanical forcing from thermally damped 
Kelvin waves. The expression for the mean flow acceleration at the equator is, 
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Here, zq is taken to be 14 km (138 mb), B is the vertical momentum flux at z$, a is the 
thermal damping rate as a function of height, c is the zonal phase speed specified as 50 
ms -1 , and k — 1 is the zonal wavenumber. We specified B to be .0105 m 2 s~ 2 which is a 
factor of 1.5 larger than that used by Gray and Pyle [1987], We also used a thermal damping 
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rate which was slightly slower than the ’fast’ rates of Dunkerton [1979], with a minimum 
e-folding damping time of 2.7 days near 50 km. We found that specifying these values of B 
and a resulted in model simulated tracer fields most similar to observations. The latitudinal 
distribution of the forcing is [Gray and Pyle, 1987], 


Fkw = 


dv>eq 

dt 


exp 


—2 $la<f> 2 
c — u 


(26) 


where a and are the radius and rotation rate of the earth, and <f) is latitude. The acceleration 
Fkw is then used in the forcing of the streamfunction in equations (5)-(6). 

Figure A.7a-f shows the wave driving and zonal mean circulation at the equator. The 
zonal mean wind shows the well known SAO signature above 10 mb, with a phase reversal 
between the upper stratosphere and upper mesosphere. Maximum Fkw occurs near the 
stratopause, with values of 2-3 ms~ l day~ l . There is a semiannual oscillation in Fkw, with 
significant values extending into the upper mesosphere during January and July, coincident 
with the weakest westerly u in the lower mesosphere. Weakest westerly accelerations in the 
upper stratosphere and lower mesosphere occur during the equinoxes, when the westerlies 
are at a maximum. This is consistent with the modeling study of Dunkerton [1979] and the 
observations of Hirota [1978]. There are also strong semiannual signals in the total heating 
rate and gravity wave forcing (Few), especially in the middle and upper mesosphere. The 
resultant w* field in the upper mesosphere shows a semiannual cycle with maximum upward 
(downward) motion during the solstices (equinoxes). However, the upper stratosphere and 
lower mesosphere exhibits a dominant annual cycle with only a weak semiannual component, 
probably due to the competing influences of the heating and wave driving. The seasonal 
variation in K zz also shows a strong semiannual signal above 80 km with maximum mixing 
during the equinoxes. This semiannual signal reverses phase in the 40-60 km region, although 
with significantly weaker modulation. 

Finally, we show the total wave driving ( Fjot = Fpw + Fqw + Fkw) for January in 
Figure A. 8. Gravity wave drag dominates the total throughout the upper mesosphere, with 
the planetary wave effects being a greater contribution to the total in the stratosphere. The 
signal in the middle to high latitude troposphere is due to the combination of planetary and 
baroclinic eddies. 


A . 2. d. Troposphe ric/ Stratospheric K zz 

We prescribe a baseline K zz which is added to the K zz computed from the gravity 
wave parameterization (section A.2.b). In the troposphere and lower stratosphere at and 
below 30 km, the background is based on the climatological temperature lapse rate, and is 
determined as follows. We specify the largest value of K zz to be 100 m 2 s for the steepest 
lapse rate observed in the lower troposphere (-8 K/km). Following the analyses of Mote et 
al. [1998] and Hall and Waugh [1997b], we specify the minimum K zz to be ,01 m 2 s~ l in 
the tropical lower stratosphere where the lapse rate is most stable (~ +4K/km). The K zz 
values at each grid point throughout the year are then interpolated from these extreme values 
using the corresponding lapse rate. This methodology assumes that larger mixing occurs in 
regions of weaker static stability, given that the latter is indicative of stronger convective 
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overturning. Here we have not accounted for the effects induced by instability in the vertical 
wind shear given the difficulty in quantifying the diffusion due to this process. Above 30 km, 
the background K zz is increased from its value at 30 km to .3 m 2 s~ l in the mesosphere. As 
seen in Figure A.5d, this K zz profile mimics the large vertical mixing in the troposphere with 
a sharp decrease in mixing across the tropopause into the lower stratosphere. 


A. 3. Numerical Advection Scheme 

We have implemented a new numerical advection algorithm replacing the previous scheme 
based on Prather [1986]. The new algorithm is a two dimensional version of the flux form 
scheme developed by Lin and Rood [1996] (hereafter LR96) currently used in the GSFC 3-D 
chemistry and transport model [Douglass et ai, 1996]. The new scheme is mass conserving 
and utilizes an upstream piecewise parabolic method (PPM) [Colella and Woodward, 1984; 
Carpenter et al., 1990]. We use the fully monotonic PPM so that no new minima or maxima 
are generated by advection. The algorithm of LR96 also includes cross terms to account 
for errors produced by operator splitting, i.e., the successive application of meridional and 
vertical advection operators. A time step of 12 hours is used for the advection of constituents. 
The PPM is highly accurate (to 4th order) and preserves sharp tracer gradients quite well, 
while exhibiting very little numerical diffusion. This was found to be especially important for 
simulations in regions of strong wind shear and sharp gradients in the tracer field, such as 
water vapor near the tropical tropopause (section 3.2). 
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Figure Captions 

Figure 1 . Zonal mean CH4 (ppbv) from the model simulation (solid line) and the 
HALOE+CLAES climatology (dashed line) for January, April, July, and October. The 
contour interval is 300 ppbv, and the 150 ppbv contour is also included. 

Figure 2. Zonal and monthly averaged month-height sections of CH4 (ppbv) from the 
model simulation and the HALOE+CLAES climatology for the northern and southern polar 
regions. Note that the plots for the NH have been shifted by 6 months to facilitate the visual 
comparison. The contour interval is 100 ppbv. 

Figure 3 . Model simulation of CH 4 in ppbv (color) for January and September, along with 
streamlines depicting the sense of the model residual circulation (v*,w*). 

Figure 4. Model simulation of CH4 as in Figure 3 for January and September, (solid contours) 
along with the model K yy fields depicted in color. The contour interval for CH4 is 100 ppbv, 
and the K yy values are in 10 8 cm 2 sec" 

Figure 5 . Zonal and monthly averaged month-latitude sections of CH4 (ppbv) from the model 
simulation and the HALOE+CLAES climatology at 2.2 mb. The contour interval is 100 ppbv. 

Figure 6 . Zonal mean H2O (ppmv) from the model simulation (solid line) and the 
HALOE+MLS climatology (dashed line) for January, March- April, July, and September- 
October. The HALOE version 18 data averaged over 1991-1997 is included above 60 km. The 
contour interval is 1 ppmv. 

Figure 7 . Zonal mean time series of HALOE H 2 0 (ppmv) averaged from 10°S and 10°N 
for 1993 to 1995. Also shown is the one year climatological model simulation of H 2 0 at the 
equator repeated for three years. Water vapor is shown in color and in the overlaid contours 
in intervals of .2 ppmv. 

Figure 8 . Zonal mean time series of the quantity 2CH 4 + H 2 0 in ppmv from HALOE 
averaged from 10°S and 10°N for 1993 to 1995. Also shown is the one year climatological 
model simulation of 2CH4 + H 2 0 at the equator repeated for three years. The values are shown 
in color and in the overlaid contours in intervals of .2 ppmv along with the 6.7 ppmv contour. 

Figure 9 . Vertical profiles of the fractional amplitude and phase lag of the annual cycle in 
2CH 4 + H 2 0 at the equator from the model (solid line) and HALOE data (dashed-asterisk 
line). These quantities were obtained by Fourier analysis of the values shown in Figure 8. The 
amplitude has been normalized to the value at the tropical tropopause (100 mb), and the 
phase lag (years) is defined to be zero at the tropical tropopause. 

Figure 10. Vertical profiles of the model dynamics in the equatorial lower stratosphere 
for January (thin solid), April (dotted), July (dashed), October (dash-dot), and the annual 
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mean (heavy solid). The dilution rate (I0~ 6 sec~ l ) and its reciprocal, the dilution time scale 
(months), are estimates of the horizontal in-mixing of mid-latitude air and are based on 
the model K yy values (see text for details). Also included are the residual vertical velocity 
(mm /sec) and the vertical diffusion (m 2 /sec). 

Figure 11. Seasonal profiles of zonal and monthly mean H 2 O from the model (solid) and 
HALOE data (dashed-asterisk) for 80 km (.01 mbar) at 45°N, equator, and 45°S. 

Figure 12. Same as Figure 11 for 65 km (.1 mbar). 

Figure 13. Zonal and monthly mean total ozone from Nimbus-7 TOMS version 7 data 
averaged over 1988 to 1992 (top), the model simulation corresponding to 1990 Cl y loading 
(middle), and the difference, model minus TOMS (bottom). The contour interval is 20 Dobson 
units (DU). 

Figure 14. Time dependent model simulations of carbon 14 every six months between 
January 1964 and July 1966. The contour interval is 50 mixing ratio units, defined as 10 5 
atoms of 14 C per gram of dry air [Kinnison et aL, 1994]. 

Figure 15, Vertical profiles of carbon 14 at 31°N every six months between January 1964 
and July 1966. Plotted are the time dependent model simulations using the new transport 
algorithm (solid), the previous 1995 transport algorithm (dotted), and observations (Kinnison 
et aL [1994] - dashed-asterisk). Values are in mixing ratio units, defined as 10 5 atoms of 14 C 
per gram of dry air. 

Figure 16. Vertical profiles of strontium 90 at 9°N for six time periods between January 
1965 and October 1966. Plotted are the time dependent model simulations using the current 
transport algorithm (solid), the previous 1995 transport algorithm (dotted), and observations 
(Kinnison et aL [1994] - dashed-asterisk). Units are proportional to mixing ratio as in Kinnison 
et aL [1994]. Model simulations include a settling velocity for strontium-90 as discussed in the 
text. 

Figure 17. Vertical profiles of strontium 90 at 34°S for six time periods between April 
1965 and October 1966. Plotted are the time dependent model simulations using the current 
transport algorithm (solid), the previous 1995 transport algorithm (dotted), and observations 
(Kinnison et aL [1994] - dashed-asterisk). Units are proportional to mixing ratio as in Kinnison 
et aL [1994]. Model simulations include a settling velocity for strontium-90 as discussed in the 
text. 

Figure 18. Mean age of air (years) derived from time dependent model simulations of SF6 
for January, April, July, and October. The age is taken relative to the global mean value at 
the surface. The contour interval is 0.5 years. 

Figure 19. Age of air as a function of latitude for October/November 1994 and 
January /February 1996. The ages are derived from SF 6 using ER-2 measurements at 19-21 
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km (triangles) and model simulations at 20 km (lines). Included are model simulations using 
the previous 1995 model transport (dotted line), the current 1998 transport (solid line), and 
the 1998 transport with a photochemical loss identical to CFC-115 imposed on SFe (dashed 
line). The age is taken relative to the global mean value at the surface. 

Figure 20. Vertical profiles of age of air from the model and balloon data taken at the 
latitudes and seasons indicated. The observations are from OMS SF6 (triangles) and the 
balloon SF6 measurements of Harnisch et al. [1996] (asterisks and squares). Included are model 
simulations using the previous 1995 model transport (dotted line), the current 1998 transport 
(solid line), and the 1998 transport with a photochemical loss identical to CFC-115 imposed 
on SFe (dashed line). The age is taken relative to the global mean value at the surface. 

Figure 21. Age of air derived from SF6 as a function of latitude for October/November 1994 
and January/February 1996. The ER-2 observations (triangles) are identical to Figure 19. 
The model simulations are for 20 km and include: the standard T998 model 1 as in Figure 19 
(solid line); the standard T998 model 1 with the minimum lower stratospheric Ii zz increased to 
.1 m 2 s~ l (dashed line); and the standard ’1998 model’ with the minimum lower stratospheric 
K zz increased to 1 m 2 s~ l (dash-dot-dot-dot line). No mesospheric loss is imposed on SF6« 

Figure 22. Vertical profiles of age of air from the model and balloon data taken at the 
latitudes and seasons indicated. The observations are identical to Figure 20 (OMS - 
triangles; balloon SF6 measurements of Harnisch et al. [1996] - asterisks and squares). The 
model simulations include: the standard T998 model’ as in Figure 20 (solid line); the standard 
’1998 model 1 with the minimum lower stratospheric K zz increased to .1 m 2 s~ l (dashed line); 
and the standard ’1998 model’ with the minimum lower stratospheric K zz increased to 1 m 2 s~ l 
(dash-dot-dot-dot line). No mesospheric loss is imposed on SFq. 

Figure 23. Latitude-height cross sections of the model annual cycle of CO 2 : (a) amplitude 
(ppmv), defined as one-half of the difference from peak to trough of the seasonal variation; (b) 
phase (month of maximum). 

Figure 24. Time series of the model CO 2 (ppmv) at 55°N for 1994 and 1995 plotted at the 
indicated pressure levels. 

Figure 25. Time series of CO 2 (ppmv) for 1992 to 1997 from the model simulation at the 
tropical tropopause (solid line), the average of the surface time series at Mauna Loa (19°N) 
and Samoa (14°S) with a two month time lag (dashed line), and averages of ER-2 CO 2 
measurements taken just above the tropopause [Boering et al, 1996] (triangles). 

Figure 26. Time series of CO 2 (ppmv) for 1994 and 1995 from ER-2 observations (top 
row) and the model (bottom row) for three latitude ranges: tropics (6°S-12°N), subtropics 
(12°S-30°N), and midlatitudes (30°S-48°N). The observations are least squares fits to the 
seasonal cycle and trend and are adapted from Strahan et al. [1998]. The data have been 
binned by N 2 O and encompass three potential temperature ranges: 370-410K (solid line), 
395-435K (dashed line), and 440-480K dashed-dotted line). See text for details. 
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Figure 27. Time series of C0 2 (ppmv) for 1994 and 1995 as in Figure 26, but plotted only 
for 370-4 10K for the three latitude regions indicated. 

Figure 28. Vertical profiles of C0 2 (ppmv) versus potential temperature (K) from the model 
(solid line) and ER-2 measurements (plus sign) averaged over 10°S to 10°N for the month/year 
indicated. Approximate altitude is shown on the right hand axes. 

Figure 29. Vertical profiles of C0 2 (ppmv) versus potential temperature (K) from the model 
(solid line) and ER-2 measurements (plus sign) averaged over 30°N to 40°N and 70°S to 60°S 
for the month/year indicated. Approximate altitude is shown on the right hand axes. 

Figure 30. As in Figure 29, but with C0 2 (ppmv) plotted versus N 2 0 (ppbv). 

Figure A.l. Latitude-height sections of zonal mean temperature (a), zonal wind (b), and 
residual vertical velocity (c) and horizontal velocity (d) from the model transport. The contour 
intervals are 10A" for temperature; 10 mj sec for zonal wind; 0.5 cm/ sec for w* including the 
contours for ±.2, ±.l, and ±.05cra/,sec; and 2 m/sec for v* including the contours for ±1,±.5, 
and ±.2 mj sec. Negative values are shaded. 

Figure A. 2. Net diabatic heating rates for January. The contour intervals are 2 Kjday and 
include the contours for ±1, and ±.5 Kjday. Negative values are shaded. 

Figure A. 3 . Contributions to the total heating rate for January from, latent heating (a), 
gravity wave heating effects seperated (b-c), and the net gravity wave heating (d). See text for 
details. Contour intervals are .5 K /day for latent heating, and 2 K/day including the contours 
for ±1, and ±.5 Kjday for the gravity wave heating. Negative values are shaded. 

Figure A. 4 . The total heating rate for January. Contour intervals are IK j day and include 
the contours for ±1, and ±.5K /day. Negative values are shaded. 

Figure A. 5 . Eliassen-Palm flux divergence due to planetary and synoptic scale waves (a) and 
gravity waves (b) for January. Also plotted are the associated model eddy diffusion coefficients 
K yy (c) and K zz (d). The contour intervals are, for (a): 2 m/ sec/day including ±1 and 
±.5 m/sec/day (shaded corresponds to values < -2 mj sec/day ); (b): 20 mj sec/day including 
±10, ±5, and ±1 m/sec/day (negative values are shaded); (c): lxl0 lo cm 2 /sec including 
2xl0 8 , lxlO 9 , and 5x10 9 cm 2 /sec (shaded corresponds to K yy > lxl0 10 cm 2 /sec); and, (d): 
2x10 ^cm 2 /sec including 200,300, 1000, lxlO 4 , 5xl0 4 , lxl0 5 cm 2 /sec. 

Figure A. 6. Month-latitude sections of the model Ii yy fields in color, along with u overlaid 
in solid contours, for 68 mbar (a), 12 mbar (b), 1.7 mbar (c), and .1 mbar (d). The contour 
interval for u is 10 m/sec, and the K yy values are in 10 8 cra 2 /sec. 

Figure A. 7. Month-height sections of the model dynamical fields at the equator. Shown are 
u (a), the total heating rate (b), the Eliassen-Palm flux divergence due to Kelvin waves (c) 
and gravity waves (d), the residual vertical velocity, w* (e), and the vertical diffusion, K zz (f). 
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The contour intervals are, for (a): 10 m/ sec ; (b): 2 K (day including ±1 and ±.5K /day; (c): .5 
m/ sec/ day including .1 m/sec/day ; (d): 2 mj secj day including ±1 and db.5 m/sec/day, (e): 
.2 cm/ sec including d=.l, ±.05, ±.02, ±.01cra/.sec; and, (f): .02, .1, .2,. 5, 1,2,5, 10,15,20, and 
30 m 2 j sec. Negative values are shaded for all plots. 

Figure A*8. The total Eliassen-Palm flux for January. Contour intervals are 20 m/sec/day 
and include the contours for ±1,±2,±5, and ±10 m/sec/day. Negative values are shaded. 
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